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ABSTRACT 


This  implementation  of  a  Legendre-Gauss-Lobatto  Pseudospectral  (LGLP)  al¬ 
gorithm  takes  advantage  of  the  MATLAB  Graphical  User  Interface  (GUI)  and  the 
Optimization  Toolbox  to  allow  an  efficient  implementation  of  a  direct  solution  tech¬ 
nique.  Direct  solutions  techniques  solve  optimal  control  problems  without  solving  for 
the  optimality  conditions.  Using  the  LGLP  method,  an  optimal  control  problem  is 
discretized  into  a  Nonlinear  Program  (NLP)  and  solved  using  an  NLP  solver,  avoid¬ 
ing  the  problems  of  deriving  the  conditions  of  optimality  and  solving  the  resulting 
boundary  value  problem.  The  MATLAB  GUI  implementation  solves  optimal  control 
problems  without  requiring  knowledge  of  the  specific  implementation  of  the  LGLP 
method.  The  GUI  completes  the  discretization  of  the  problem  and  solves  the  resulting 
NLP  using  a  Sequential  Quadratic  Programming  Algorithm.  The  GUI  will  convert 
any  optimal  control  problem  with  fixed,  free  or  optimal  final  time,  a  Mayer,  Lagrange 
or  Bolza  cost  function,  constrained  or  unconstrained  controls,  with  or  without  state 
inequalities,  and  point  inequalities  into  a  parameter  optimization  problem  and  re¬ 
turns  a  solution.  The  GUI  creates  a  function  file,  output  file,  binary  save  file,  and 
optimization  script  to  allow  full  access  to  the  strength  of  the  LGLP  method  from  the 
GUI  or  the  command  line.  No  prior  knowledge  of  the  LGLP  algorithm  is  assumed  or 
necessary. 
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I. 


INTRODUCTION 


Optimal  Control  Theory  is  one  of  the  modern  applications  of  the  Calculus 
of  Variations.  Motivated  by  the  development  of  modern  computers  and  the  space 
program  during  the  1960s,  the  theory  of  the  calculus  of  variations  has  been  used 
to  create  methods  for  designing  modern  systems.  In  1962  Pontryagin’s  minimum 
principle  gave  the  theoretical  basis  for  solving  optimal  control  problems.  [Ref.  2] 

The  calculus  of  variations  was  developed  to  describe  situations  that  nature  had 
already  optimized.  The  theory  seeks  to  minimize  or  maximize  the  value  of  integrals. 
A  classic  example,  known  as  the  brachistochrone  problem  [Ref.  1],  is  to  minimize 
the  time  a  bead  takes  to  slide  down  a  wire  between  two  points.  In  this  problem, 
the  shape  of  the  wire  can  be  directly  changed  to  decrease  the  time  the  bead  requires 
to  move  along  the  wire.  In  this  type  of  problem,  the  physical  states  can  be  directly 
altered.  The  system  can  be  represented  as  finding  y,  the  shape  of  the  wire,  so  that  the 
time,  expressed  as  an  integral  of  the  physical  states,  can  be  minimized.  The  following 
expression 

^(x,  y)=  /  f(x,y,y)dx  (1.1) 

J  a 

represents  the  cost  function  to  be  minimized.1 

Optimal  Control  Theory  is  a  generalization  of  the  calculus  of  variations,  but 
the  control  problems  have  unique  characteristics  that  identify  this  class  of  problem. 
Optimal  Control  problems  select  from  available  controls  to  optimize  the  performance 
of  a  dynamic  system.  These  problems  focus  on  altering  forces  at  our  disposal  to 
control  a  process  and  optimize  some  performance  measure  of  the  system  [Ref.  1]. 
The  controls  can  be  altered,  but  the  state  of  the  system  may  not  be  directly  altered. 
If  y  =  u  is  used  as  a  state  equation  and  substituted  into  equation  (1.1),  the  calculus 

throughout  this  thesis  boldface  notation  will  be  used  to  distinguish  n  x  1  vectors  or  functions 
of  n  x  1  vectors. 
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of  variations  problem  can  be  converted  into  the  following  optimal  control  problem, 

rb 

J(x,y,u)  =  j'  f(x,y,u )dx.  [Ref.  1]  (1.2) 

The  emphasis  in  optimal  control  problems  is  on  the  controlling  process.  The  ad¬ 
justable  variables  are  the  control  variables,  not  the  state  variables. 

Optimal  control  problems  also  add  constraints  to  the  calculus  of  variations. 
Controls  might  have  physical  limitations,  such  as  the  force  from  an  aircraft  engine 
or  the  length  of  a  robotic  arm,  and  these  constraints  are  modelled  into  the  dynamic 
system  that  must  be  optimized.  These  types  of  problems  in  applied  mathematics 
lead  to  engineering  applications  where  controlling  devices  can  be  used  to  guide  the 
state  equations  in  an  optimal  way.  The  historical  motivations  for  the  field  were  the 
space  and  rocket  programs  apd  problems,  such  as  the  lunar  landings.  As  the  field 
developed,  it  has  been  applied  to  financial,  robotic,  artificial  intelligence  and  other 
engineering  applications. 

Optimal  Control  Theory  seeks  to  control  a  system  while  minimizing  or  maxi¬ 
mizing  a  cost  function.  The  cost  function  will  reflect  the  change  in  some  state  of  the 
system,  some  measurable  amount  of  time,  or  the  amount  of  control  effort  required 
for  the  system  to  reach  a  desired  end  state.  Finding  optimal  functions  is  the  goal  of 
optimal  control  theory.  This  differs  from  standard  parameter  optimization  because 
with  optimal  control  theory,  optimization  is  performed  over  a  continuous  function 
space. 

In  optimal  control  problems,  dynamic  constraints  are  added  to  the  cost  func¬ 
tion  to  be  optimized.  The  state  equations  and  inequality  constraints  on  the  controls 
and  states  are  adjoined  to  the  cost  function  through  the  use  of  Lagrange  multipliers. 
The  Hamiltonian  is  used  to  simplify  the  equations  used  in  this  minimization.  The 
Hamiltonian  represents  the  sum  of  potential  and  kinetic  energies  used  in  a  dynamic 
system.  The  first  variations  of  the  cost  function  components  are  set  to  zero  to  ensure 
a  local  minimization  of  the  cost  function.  The  equations  that  result  are  the  necessary 
conditions  for  optimality  and  are  referred  to  as  the  Euler-Lagrange  equations.  An 
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example  of  an  augmented  cost  function,  Ja ,  is 

3 A  =  ;  {£(x,u,<)  +  Ar  [f(x,  u,<)  —  x]  +  nTg(u,t)\  dt,  (1.3) 

with  a  dynamic  constraint  of  x  =  f(x,  u,t),  initial  condition  of  x(t0)  =  x0,  and 
equality  constraints  of  g(u,  t)  —  0.  This  augmented  Bolza  cost  function  has  a  final 
time  cost,  and  an  integral  cost,  ///  £(x,  u,  t).  which  are  adjoined  by  the 

state  constraints.  The  variables  A  and  fi  are  Lagrange  multipliers.  The  augmented 
Hamiltonian  is  defined  as 

U  =  C  +  XTf  +  tiTg.  (1.4) 

The  optimality  conditions  are2 


The  solutions  to  the  derived  optimality  conditions  often  rely  on  numerical  methods 
for  solving  two  point  boundary  value  problems. 

Optimal  control  algorithms  can  be  divided  into  two  types,  direct  and  indirect 
methods.  Most  texts  on  optimal  control  theory,  [Ref.  1],  [Ref.  2],  [Ref.  6],  [Ref.  14], 
focus  on  indirect  numerical  techniques.  Indirect  methods  rely  on  the  necessary  op¬ 
timality  conditions  derived  from  the  minimum  principle.  The  necessary  conditions 
are  developed  and  then  solved  to  find  an  optimal  trajectory.  The  solution  normally 

2The  gradient  and  Jacobian  used  to  develop  the  optimal  conditions  will  be  defined  in  the  next 
section. 


3 


requires  solving  a  nonlinear  two  point  boundary  value  problem  (BVP),  and  solving  for 
the  Lagrange  multipliers  and  costates.  Nonlinear  BVPs  normally  do  not  have  a  closed 
form  solution  and  theorems  do  not  exist  that  guarantee  existence  and  uniqueness  for 
all  BVPs.  This  lack  of  an  analytic  solution  requires  the  use  of  numerical  methods 
to  solve  the  BVP.  The  most  common  indirect  method  in  use  today  is  an  indirect 
shooting  algorithm  [Ref.  3],  Direct  methods  discretize  the  continuous  problem  into 
a  parameter  optimization  problem  and  then  solve  the  resulting  nonlinear  optimiza¬ 
tion  problem.  Direct  methods  do  not  explicitly  employ  the  necessary  conditions  for 
optimality,  but  the  results  can  be  checked  using  the  optimality  conditions  from  the 
calculus  of  variations. 

An  advantage  of  indirect  methods  is  that  they  directly  solve  the  adjoint  dif¬ 
ferential  equations  for  the  conditions  derived  from  the  minimum  principle  and  the 
transversality  conditions  to  ensure  optimality.  The  first  derivative  of  the  Hamilto¬ 
nian  is  set  to  zero  and  the  resulting  system  of  ordinary  differential  equations  is  solved. 
This  advantage  leads  to  several  nontrivial  problems: 

1.  Necessary  conditions  must  be  derived  analytically. 

2.  Region  of  convergence  for  root  finding  algorithms  required  to  solve  the  BVP 
may  be  very  small  requiring  a  very  “good”  initial  guess. 

3.  Path  inequalities  may  require  solving  for  constrained  and  unconstrained  sub 
arcs  before  each  iteration  begins. 

4.  When  setting  the  gradient  to  zero,  analytic  expressions  for  the  gradient  must 
be  calculated. 

There  are  several  indirect  numeric  methods  for  optimal  control  problems.  One 
such  method  is  the  neighboring  extremal  method  or  shooting  method  [Ref.  6].  Starting 
values  for  the  state  vector,  a?(t0),  and  adjoint  vector  A(t0),  are  guessed  and  then  the 
state  equations 

x(t)  =f(x,u,t)  (1.8) 
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and  costate  equations 


(1.9) 


dx 

are  integrated  from  to  to  tf  with  the  control  history 

dH(x,u,  \,t) 

du 

The  optimal  trajectory  that  results  from  the  solution  of  the  differential  algebraic 
equation  does  not  necessarily  satisfy  the  end  conditions.  The  starting  values  of  the 
adjoint  variable,  A,  can  be  perturbed  in  various  ways  to  meet  the  end  conditions. 
Generally,  shooting  methods  are  very  sensitive  to  variations  in  the  initial  guess. 

A  second  indirect  method  is  a  gradient  method,  the  method  of  steepest  descent. 
With  gradient  methods  the  dynamic  system  equations  are  solved  exactly  at  each 
iteration  and  the  control  is  slightly  perturbed.  The  control  history  is  adjusted  at 
each  step  to  further  reduce  the  cost.  The  state  equation  (1.8) 

x(t)  =  f(x,u,t)  (1.11) 

is  integrated  from  t0  to  tf  to  obtain  x(t)  for  a  guessed  u(i),  t  €  \to,tf\.  The  cost 
sensitivity  matrices  and  Lagrangian  gradients  are  then  evaluated.  The  adjoint  vector 
is  integrated  backwards  to  obtain  A (f),  thus  determining  Hu.  7iu  may  not  start  close 
to  zero,  but  approaches  zero  with  each  iteration.  This  method  of  convergence  makes 
stopping  criteria  very  important  to  this  method. 

Direct  methods  have  been  developed  to  avoid  solving  the  optimality  condi¬ 
tions.  Instead  of  starting  by  solving  for  the  optimality  conditions,  direct  methods 
begin  by  discretizing  the  problem.  The  states,  controls  or  both  are  discretized  and 
the  continuous  problem  changes  into  a  discrete  problem  before  the  cost  function  is 
optimized.  This  reduces  the  optimal  control  problem  to  a  parameter  optimization 
problem  that  can  be  solved  by  Nonlinear  Programming  (NLP)  solvers.  The  accuracy 
of  the  solution  depends  on  the  discretization  producing  a  problem  that  accurately 
represents  the  continuous  problem.  With  the  current  advances  in  NLP  solvers,  a 


=  0. 


(1.10) 
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properly  discretized  problem  will  result  in  solutions  that  are  increasingly  more  accu¬ 
rate.  In  the  1960s,  Newton’s  method  was  used  to  solve  direct  problems  and  the  size  of 
the  problems  was  limited  to  n  =  m  =  10,  with  n  states  and  m  controls.  In  the  1970s, 
the  use  of  quasi-Newton  methods  and  reduced  gradient  methods  increased  the  size 
to  n  =  m  <  100.  Current  advances  in  computing  with  sparse  matrices  in  Numerical 
Linear  Algebra  have  allowed  solving  problems  with  n  =  m  =  10,000  [Ref.  3],  These 
advances  have  resulted  in  continued  interest  and  research  in  direct  methods. 

This  thesis  focuses  on  a  direct  method  of  solving  optimal  control  problems. 
It  relies  on  a  Legendre-Gauss-Lobatto  Pseudospectral  (LGLP)  algorithm  which  was 
first  utilized  in  [Ref.  4]  for  solving  problems  in  optimal  control  theory.  This  algorithm 
discretizes  the  problem  and  converts  the  problem  into  an  NLP.  Further  developments 
by  Professor  Fahroo  and  Professor  Ross  [Ref.  5],  implemented  the  algorithm  and 
performed  costate  estimation  in  MATLAB. 

The  original  research  in  this  project  involved  creating  a  Graphical  User  Inter¬ 
face  (GUI)  in  MATLAB  that  allows  optimal  control  problems  to  be  solved  using  the 
LGLP  method  without  requiring  a  thorough  understanding  of  the  method.  The  end 
state  of  this  project  is  a  GUI  that  only  requires  understanding  of  the  basic  optimal 
control  problem  formulation  to  allow  solution  of  a  wide  variety  of  problems.  The 
LGLP  code  relies  on  codes  developed  by  Professor  William  Gragg  to  solve  Numerical 
Linear  Algebra  problems  associated  with  the  algorithm.  In  addition  to  a  knowledge 
of  MATLAB  and  Numerical  Linear  Algebra,  a  basic  understanding  of  the  Calculus 
of  Variations  and  Optimal  Control  Theory  is  required  for  this  thesis. 
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II. 


OPTIMAL  CONTROL  PROBLEMS 


A  dynamic  system  that  represents  an  optimal  control  problem  is  mathemat¬ 
ically  described  by  a  system  of  ordinary  differential  equations.  The  optimal  control 
problem  is  defined  by  state  equations,  boundary  conditions  on  the  state  variables, 
equality  constraints  and/or  inequality  constraints  on  the  state  variables,  constraints 
on  the  controls,  and  a  cost  function  to  be  optimized. 

Definition:  Let  s( y)  be  a  scalar  function  of  y  =  [yi . . .  yn]T ■  The  gradient  of  s 
with  respect  to  y  is  defined  as  the  n  x  1  vector: 


ds  _ 

§7(y)  = 


t(y) 

fe(y) 


— (y) 


L  dyn 


(2.1) 


Definition:  If  a(y)  is  an  n  x  1  vector  function  of  y  (an  m  x  1  vector),  the 
Jacobian  matrix  is  defined  as  the  n  x  m  matrix: 


>> 

1 _ 

fe(y)  - 

fc(y) 

d[a(y)]  _ 

&(y) 

&(y)  - 

fe(y) 

dy 

• 

l 

• 

.fe(y) 

fe(y)  ••• 

fe(y) 

(2.2) 


A.  COST  FUNCTIONS 

The  cost  function  to  be  optimized  is  a  function  of  the  states,  controls  and  final 
time.  There  are  three  types  of  cost  functions:  Mayer,  Lagrange  and  Bolza.  A  Mayer 
cost  function  is  of  the  form: 


J(u,x.,tf)  =  M(x,tf).  (2.3) 

The  Lagrange  type  of  cost  function  has  an  integral  cost.  It  has  the  form: 

k7(u,x,</)=  f  '  £(x,u)dt.  (2.4) 

Jto 


7 


The  Bolza  cost  function  incorporates  both  the  Mayer  and  Lagrange  forms: 

J (u,  x,tf)  =  M (x,  tf)+  f  !  £(x,  u)df .  (2.5) 

It  can  be  shown  that  the  three  types  of  cost  functions  are  equivalent  and  any  type  of 
cost  function  may  be  transformed  into  either  of  the  other  types  [Ref.  6]. 

B.  STATE  EQUATIONS  AND  CONSTRAINTS 

The  solution  of  the  problem  is  found  by  determining  the  control  function  u(t), 
and  the  corresponding  state  trajectory  x(f),  that  minimizes  the  cost  function.  The 
problem  has  a  state  vector  x  €  Rn  and  a  control  vector  u  <E  Rm.  The  dynamic 
equations  for  the  states  are  expressed  as: 

*(*)  =  f(x(*),u(i)),  t  e  [t0,tf]  (2.6) 

with  initial  and  final  boundary  conditions: 

lMx(*o),<o]  =  0,  (2.7) 

^/[x(*/M/]  =  0,  (2.8) 

where  E  ft?  with  p  <  n  and  ^ /  €  Rq  with  q  <n.  The  problem  may  have  inequality 
constraints  on  the  controls  which  are  formulated  as 

«[«(*)]  <  0  (2.9) 

where  g  €  RT  and  |g  has  full  rank. 

C.  THE  ADJOINT  DIFFERENTIAL  EQUATIONS 

Once  the  problem  is  mathematically  formulated,  the  principles  from  calculus 
of  variations  are  used  to  derive  necessary  and  sufficient  conditions  for  optimality.  We 
adjoin  the  state  equations  and  constraints  to  the  cost  function  to  create  an  augmented 
cost  function.  For  simplicity  in  the  presentation  of  the  adjoint  differential  equations 
and  the  following  optimality  conditions,  we  make  the  simplifying  assumption  that 

sN*)]  =  o.  (2.10) 
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For  a  complete  treatment  of  the  inequality  case,  reference  [Ref.  14].  The  augmented 
cost  function  is 

Ja  =  v7  +  ATf  +  tiTg  +  v*ip0  +  tfipf.  (2.11) 

Using  Lagrange  multiplier  theory,  \{t)  €  Rn,  fi{t)  €  Rr,  u0 (t)  6  RP,  with  p  <  n, 
and  *//(<)  €  J?9,  with  q  <  n,  are  the  Lagrange  multipliers.  The  variable  f  repre¬ 
sents  the  state  equations,  ip 0  represents  initial  time  constraints,  g  represents  control 
constraints,  and  ip j  represents  final  time  constraints.  Since  the  state  equation  is 

x  =  f(x,u),  (2.12) 

we  adjoin  the  constraint 

f(x,u)-x  =  0  (2.13) 

to  the  cost  function  resulting  in 

Ja=J  +  AT(f  (x,  u)  -  x)  +  fiTg(u)  +  i'oV’o  +  J'jqfr  /  (2.14) 

Using  a  Bolza  cost  function,  the  augmented  cost  function  is 

Ja  =  M(x.tj)  +  vliia  +  +  J^'  [£(x,u)  +  Ar(f(x,u)  -x)  +  /jTg(u)]  it. 

(2.15) 

Using  the  Hamiltonian  in  the  cost  function  will  simplify  the  equations.  The  Hamil¬ 
tonian  is  defined  as 


■H(x,  A,  u,  ft)  =  £(x,  u)  +  Arf(x,  u)  +  fiTg(u). 


(2.16) 


Thus  by  expanding  the  cost  function, 

r*f 


Ja  -  M(x.,tf)+ulip0+u'jtpJ+J^C(x,u)dt+J^ffiT^dt+ J  *  Arf(x,u)  dt- J  f  XTxdt, 

‘(2.17) 

and  substituting  the  Hamiltonian,  we  have 


Ja  =  M(x,</)  +  VqiPq  +  vT}ip}  +  £  f  U (x,  A,  u,  fi)  dt  -  jf  '  \Txdt.  (2.18) 
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For  simplicity,  the  augmented  Bolza  cost  function  is  expressed  as 


Ja  =  M[-]  +  o  +  i/JV,  +  jf  r  W[*]  -  j*f  \Txdt.  (2.19) 

Using  integration  by  parts  in  the  second  integral  we  have 

£  XTiii  =  AH'  -  £  xTxdt-  (2.20) 

This  allows  simplification  of  the  cost  function  to 

J  =  M[]  +  u^0  +  -Arx|^  +  ATx|io  +  j {*'  H[-]dt  +  j f*f  ATx  d*.  (2.21) 

D.  NECESSARY  OPTIMALITY  CONDITIONS 

The  minimum  of  the  functional,  J ,  for  the  optimal  control  u*  occurs  when 

J (u)  -  J( u*)  =  AJ  >  0  (2.22) 

for  all  controls  u  close  to  u*.  Let  u  =  u*  +  du  and 

Ai7(u*,du)  =  <L7(u*,<fa)  +  higher  order  terms.  (2.23) 

8J  is  linear  in  <Ju  and  the  higher-order  terms  approach  zero  as  the  norm  of  du 
approaches  zero.  The  increment  of  J  is 

AJ(u*,Su)  =  —(x(tf)ttf)+  *'/“<*(*/)  8xf 

+  uj  +  Stf 

rtf  ( r .  dH  1  ^ 

+  k  ( [A^  +  &T(x(^ u^’ A^’ 

r  Q'ft  1  t 

+  [du(X(t)’U(<),A(^),#<(*)’<)  (2.24) 

d'H  1 r  ) 

+  [^(*W>  “«>  A(<),  rit),  t)  6p(t)  |  dt 

+  higher  order  terms. 


By  setting  the  coefficients  of  the  differentials  equal  to  zero,  various  state  and  control 
constraints  are  obtained.1  If  the  state  equations  are  satisfied,  and  the  constraints  on 


control  are  satisfied 


then  A(f)  is  selected  so  that 


(x(t),u(t),A(i), /*(<)»<)  =  0, 


Mt)  +  fo-(x(*)>  u(*)>  *(*)>/*(*)>  *)  =  °> 


(2.25) 


(2.26) 


and  the  final  time  conditions  are 


"  T  T 

Sxf  (2.27) 

+  ^(x(</).u(</),A(</),/i(t/),</)+  "S  +  ^(x(</)></)  Stf  =  0. 

For  final  state  specified  and  final  time  free,  Sxf  =  0  and  8tj  is  free.  For  final  state 
free  and  final  time  specified,  8xj  is  free  and  8tf  =  0.  For  more  complicated  examples 
refer  to  [Ref.  2]. 

The  increment  of  J  simplifies  to 

f  dH  1^  1 

Aj7”(u*,^u)  =  Jt  <  -^-(x(t),  u(<),  A(t),  fi(t),t)  £u(f)  >  dt  +  higher  order  terms. 

(2.28) 

This  integral  is  the  first  order  approximation  of  the  change  in  "H  due  to  the  control  u, 


0^(*(O,u(<),A(*)» #*(<)> *)J  ~  ^(x(f),u*(<)  +  <Ju,A(t),^(t),«) 

-n(x(t),  u*(f),  A  (f),  t).  (2.29) 

Simplifying,  the  increment  of  J  is 

AJ(u*,Su)  =  f  f  +  8\i,\(t),n(t),t)  (2.30) 

JtQ 

—  [H(x(<),u*(f),  A(f),/x(t),f)]}  dt  +  higher  order  terms. 


1In  our  notation  the  gradient  of  a  scalar  is  a  n  x  1  vector  and  the  Jacobian  of  an  n  x  1  vector  of 
m  variables  has  dimension  n  x  m.  The  variational  Sx  is  also  n  x  1. 
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If  u*  +  Su  is  in  a  suitably  small  neighborhood  of  u*,  (||5u||  <  /?),  the  higher  order 
terms  may  be  neglected  and  the  necessary  condition  for  optimality  is 


ft(x(t),  u *(t)  +  5u,  t )  -  U{x{t),  u*(t),  A (<),  MO,  t)  >  0  (2.31) 

or 

u .*(<)  +  <Ju,  A(f ),  MO,  0  >  ^(x(t),  u*(0,  A(0,  n(t),  t).  (2.32) 

The  necessary  condition 

U(x(t), u *(f),  A(<),  MO,  0  <  ^(x(0,  u(0,  A(0,  MO,  0  (2.33) 

is  Pontryagin’s  minimum  principle.  The  minimal  principle  states  that  an  optimal 
control  must  minimize  the  Hamiltonian.  [Ref.  2] 

In  summary,  the  optimal  control  problem  is  formulated  as  follows:  find  con¬ 
trols,  u*  €  U,  the  set  of  all  possible  controls,  that  cause  the  system  described  by  Equa¬ 
tions  (2.6  -2.8)  and  (4.7)  to  follow  a  feasible  trajectory  to  minimize  the  cost  function, 
Equation  (2.5).  The  necessary  conditions  for  optimality  are  Equations  (2.25-2.28) 
and  Pontryagin’s  minimum  principle  Equation  (2.33). 

E.  NUMERICAL  SOLUTIONS 

Once  the  necessary  conditions  have  been  derived,  an  indirect  numerical  method 
may  be  used  to  solve  the  resulting  system.  The  derivation  of  the  problem  is  essential 
to  the  understanding  of  the  problem,  but  is  only  required  in  the  indirect  solution 
techniques.  Solving  the  nonlinear  BVP  problem  that  is  described  by  this  system  of 
equations  can  be  very  difficult.  To  avoid  this  computational  difficulty,  we  will  focus 
on  direct  methods.  The  use  of  a  direct  method  allows  the  solution  of  the  optimal 
control  problem  by  use  of  existing  NLP  solvers  and  does  not  require  solving  for  the 
optimality  conditions. 
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Ill 


DIRECT  METHODS 


Direct  methods  are  more  recent  developments  when  compared  with  their  in¬ 
direct  method  counterparts.  All  direct  methods  involve  deciding  upon  a  finite  set  of 
variables  and  then  iterating  by  Newton’s  method  or  another  zero  finding  algorithm  to 
numerically  solve  the  resulting  problem.  The  essential  difference  between  direct  and 
indirect  methods  lies  in  the  direct  methods  converting  the  optimal  control  problem 
into  an  NLP  and  directly  solving  the  NLP,  instead  of  deriving  the  optimality  condi¬ 
tions  and  solving  the  resulting  BVP.  The  two  most  commonly  used  direct  methods 
are  direct  shooting  and  direct  transcription  algorithms  [Ref.  3].  Differential  inclusion, 
a  technique  for  reducing  the  size  of  the  NLP  is  also  included  in  this  thesis. 

In  direct  methods,  the  Karush-Kuhn-Tucker  conditions  from  NLP  theory  are 
solved  during  the  optimization  instead  of  solving  for  the  optimal  control  necessary 
conditions.  It  is  not  surprising  that  when  the  discrete  conditions  are  taken  to  the 
limit,  determining  the  NLP  active  set  is  equivalent  to  finding  the  constrained  subarcs 
and  junction  points  in  the  continuous  optimal  control  problem. 

The  first  step  in  a  direct  method  is  to  select  a  method  to  discretize  the  problem 
and  convert  it  into  a  parameter  optimization  problem.  The  problem  is  discretized  by 
dividing  the  time  interval  into  subintervals  whose  endpoints  are  called  nodes.  These 
nodes  will  define  the  variables  for  the  NLP  problem.  The  unknowns  are  the  values 
of  the  controls,  states,  and  problem  parameters  at  the  nodes.  These  variables  are 
the  state  and  control  parameters  for  the  new  parameter  optimization  problem.  When 
the  state  equations  and  cost  function  are  expressed  in  terms  of  these  parameters,  the 
optimal  control  problem  is  reduced  to  an  NLP.  [Ref.  7] 

The  second  selection  that  defines  a  direct  method  is  the  choice  of  an  inter¬ 
polation  scheme  to  interpolate  the  values  or  time  histories  of  the  control  and  state 
variables  between  the  discrete  nodes.  Either  an  explicit  or  implicit  scheme  will  be 
chosen  to  describe  the  equations  of  motion.  An  explicit  scheme  solves  the  ordinary 
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differential  equations  (ODE)  to  develop  the  time  histories  between  the  nodes.  An 
implicit  scheme  uses  integration  rules  to  create  formulas  for  the  integrals  and  then 
adds  these  equations  as  constraints  to  the  NLP.  Many  implicit  interpolation  schemes 
use  orthogonal  polynomials  to  interpolate  the  state  and  control  values  between  the 
nodes.  One  advantage  of  using  orthogonal  polynomials  is  their  close  relationship  to 
Gauss  integration  rules,  which  results  in  highly  accurate  quadrature  rules. 

For  explicit  integration  the  value  of  x,  the  state  variable,  is  needed  to  evaluate 
the  value  of  the  function  at  each  discretization  point.  Explicit  integration  allows 
integration  of  the  state  equations  from  initial  to  final  time  in  one  pass.  With  implicit 
integration,  the  values  of  x  are  not  known  in  advance  and  as  a  result  a  predictor- 
corrector  approach  must  be  used. 

It  is  desirable  for  all  integration  rules  to  be  of  the  highest  order  possible. 
The  higher  the  order  of  an  integration  scheme,  the  smaller  the  error  from  a  single 
integration  step  taken  with  a  step  size  8.  For  example,  with  a  step  size  of  S  =  0.1,  and 
using  Euler’s  quadrature  rule  whose  local  error  is  0(82),  the  error  would  be  «  0.01. 
For  a  higher  order  rule,  such  as  Simpson’s  rule  with  a  local  error  of  0(85),  the  error 
would  be  «  0.00001  for  the  same  step  size.  This  error  dependence  on  step  size  and 
integration  rule  is  shared  by  all  the  direct  methods. 

Inequality  constraints  make  optimal  control  problems  increasingly  difficult  and 
efficiency  of  handling  such  constraints  is  very  important.  Constrained  arcs  are  often 
not  known  a  priori  and  the  junction  points  from  constrained  to  unconstrained  arcs 
may  also  not  be  known.  These  discontinuities  of  the  control  and  adjoint  variables  at 
the  junction  points  are  essentially  boundary  conditions  and  transform  the  problem 
into  a  multi-point  BVP. 

Table  I  illustrates  a  basic  algorithm  for  a  direct  method  to  convert  an  optimal 
control  problem  into  a  parameter  optimization  problem,  equivalent  to  an  NLP  [Ref.  7]. 
The  direct  methods  mentioned  here,  direct  shooting  and  direct  collocation,  both  dif¬ 
fer  in  their  choice  of  variables  for  the  NLP.  The  direct  shooting  method  discretizes 
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1 .  Divide  the  time  interval. 

2.  Choose  variables  to  be  calculated  by  interpolation. 

3.  Integrate  the  state  equations  explicitly  or  implicitly. 

4.  Solve  the  NLP. 

Table  I.  Direct  Method  Algorithm 

the  control  history,  and  direct  collocation  discretizes  both  states  and  controls.  This 
difference  tends  to  make  a  direct  collocation  NLP  larger  than  the  resulting  systems 
from  the  other  methods.  The  following  summaries  of  direct  methods  are  taken  from 
John  Betts  “Survey  of  Numerical  Methods  for  Trajectory  Optimization”  [Ref.  3]. 

A.  DIRECT  SHOOTING 

In  a  direct  shooting  method  the  variables  for  the  NLP  are  chosen  from  the 
initial  conditions,  final  conditions  and  problem  parameters.  Phases  are  chosen  at 
points  where  the  problem  is  discontinuous  and  at  boundaries  between  constrained 
and  unconstrained  arcs.  The  control  history  will  be  represented  by  a  finite  set  of 
parameters.  For  each  phase,  k ,  the  NLP  variables  are  defined  as 

X{k)  =  {x(t0),  x(</),  t0,  tf,  p}  (3.1) 

where  the  x  represents  the  state  equations  and  p  represents  a  set  of  parameters. 
[Ref.  7] 

The  control  u  can  be  explicitly  or  implicitly  represented  by  p  as  in  the  following 
examples 

u  =  pi  +  p2t  (3.2) 

or  implicitly  as 

0  =  Piu(t)  +  sin[p2«(0]-  (3-3) 
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The  total  set  of  NLP  variables  will  be 


(3.4) 

Solving  for  the  constraints  at  the  end  points  of  each  phase  provides  constraints 

g(x)  =  {^(1)(x(t0),  t0,  p),  ^(xfa),  tu  p), . . . ,  tf,  p)}  (3.5) 

where  ip^  denotes  a  function  that  describes  the  boundary  conditions.  [Ref.  3] 

The  Program  to  Optimize  Simulated  Trajectories  (POST)  developed  by  Mar¬ 
tin  Marietta  and  the  Generalized  Trajectory  Simulation  (GST)  program  by  the  Aerospace 
Corporation  are  both  implementations  of  this  method.  The  flight  avionics  on  the 
space  shuttle  also  incorporate  a  direct  shooting  algorithm  for  steering.  [Ref.  3] 

The  direct  shooting  algorithm  has  advantages  for  launch  vehicle  and  orbit 
transfer  problems.  These  types  of  problems  result  in  an  NLP  represented  by  a  small 
number  of  variables.  The  direct  shooting  algorithm  with  a  large  number  of  variables 
tends  to  lose  accuracy  due  to  propagation  of  early  errors  throughout  the  procedure. 
Another  disadvantage  is  the  cost  of  computing  the  gradient  functions.  The  system 
equations  must  be  integrated  forward  numerically  in  direct  shooting  methods.  Nor¬ 
mally,  finite  difference  codes  are  used  to  solve  for  the  gradients.  In  addition  to  the 
numeric  complexity  involved  with  solving  for  the  gradient  of  a  large  system,  there  are 
also  problems  with  the  accuracy  of  the  gradient  information.  Using  forward  difference 
methods  to  calculate  the  gradients  involves  an  error  of  the  order  of  the  step  size,  0(8). 
Using  central  difference  methods  to  calculate  the  gradients  is  twice  as  expensive  as 
forward  difference  methods,  but  the  error  is  0(82). 

A  direct  shooting  algorithm  is  contained  in  Table  II.  Within  Table  II,  x 
represents  the  state  equations,  u  represents  the  control,  and  ip  describes  the  boundary 
conditions. 

Direct  shooting  methods  are  very  efficient  when  a  detailed  mathematical  model 
is  not  required.  In  many  modern  problems,  such  as  problems  using  current  rockets 
which  do  not  have  variable  rate  thrusters  and  have  a  very  high  weight  to  foj-ce  ratio, 
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1.  Divide  the  problem  into  k  number  of  phases  based  upon  arc  constraints  and 
inequalities. 

2.  Generate  the  control  u  implicitly  or  explicitly  from  a  finite  set  of  parameters 
P- 

3.  Compute  initial  conditions,  ipk[xk(t0),pk,t0],  for  each  phase  numerically. 

4.  Solve  initial  value  problem  for  the  subarcs. 

5.  Compute  final  conditions,  tj)k[xk(tf),pk,tf],  for  each  phase  numerically. 

6.  Solve  the  resulting  NLP. 

Table  II.  Direct  Shooting  Algorithm 

detailed  mathematical  models  are  unnecessary  [Ref.  3].  The  small  number  of  variables 
in  the  NLPs  developed  by  the  direct  shooting  method  speed  up  convergence  to  a 
solution.  When  a' large  time  interval  is  used,  the  conversion  from  a  single  integration 
step  to  multiple  integration  steps  to  cover  the  time  interval  converts  the  standard 
direct  shooting  method  to  a  multiple  direct  shooting  method. 

B.  DIRECT  COLLOCATION 

A  direct  collocation  method  discretizes  both  the  control  and  the  state  equa¬ 
tions.  The  values  will  be  known  exactly  at  the  nodes  of  the  discrete  time  intervals. 
The  integration  is  completed  by  defining  the  residuals  for  each  integration  step  and 
then  driving  the  residuals  to  zero  during  the  solution  of  the  NLP  [Ref.  7].  The  NLP 
variables  are  exactly  the  values  of  the  state  and  control  equations  evaluated  at  the 
collocation  points.  When  low  order  rules  are  used,  the  method  is  called  transcription. 
When  higher  order  rules  are  used,  the  method  is  termed  collocation.  Euler  and  trape¬ 
zoidal  rules  are  used  with  transcription  methods  while  Simpson’s  rule,  Gauss- Lobatto, 
and  other  higher  order  rules  are  used  in  collocation  methods. 

A  direct  collocation  method’s  strength  lies  in  not  requiring  prior  knowledge 
or  specification  of  the  arc  sequence  for  path  inequalities.  Singular  arcs  arise  when 
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the  control  derived  is  not  uniquely  defined  by  the  optimal  conditions.  Whenever 
the  second  partial  of  the  Hamiltonian  with  respect  to  the  control  is  equal  to  zero, 
'Huu  0,  the  optimal  solution  computed  will  not  be  unique. 

The  direct  collocation  method  results  in  very  large  NLPs,  but  the  Jacobian 
and  the  Hessian  matrices  that  result  are  sparse,  allowing  for  efficient  codes  from  recent 
advancements  in  numerical  linear  algebra  to  exploit  these  matrix  structures.  It  is  not 
uncommon  for  the  matrices  involved  in  the  solution  of  a  direct  collocation  problem 
to  have  only  one  percent  non-zero  entries. 

In  addition  to  the  efficient  codes  used  with  collocation  methods,  the  high  order 
quadrature  rules  used  in  the  implicit  integration  of  the  state  equations  allow  for  a 
larger  step  size  with  an  equally  high  degree  of  accuracy  when  compared  with  the  other 
methods  listed  here.  Using  Simpson’s  rule  the  local  error  is  0(85)  and  with  Gauss 
Lobatto  the  local  error  is  0(88). 

The  implicit  integration  rules  are  added  to  the  NLP  as  nonlinear  constraint 
equations.  The  piecewise  smooth  interpolating  polynomials  used  with  direct  colloca¬ 
tion  will  satisfy  the  ODEs  exactly  at  the  collocation  points  of  each  interval.  With 
collocation,  the  states  are  unknown  and  therefore  the  equations  of  motion  must  be 
integrated  implicitly  at  each  node,  with  each  node  representing  one  integration  step. 
The  size  of  the  step  dictates  the  number  of  integration  steps  required.  A  balance  is 
desired  between  a  small  step  size  and  increasing  the  number  of  integrations.  A  direct 
collocation  algorithm  is  shown  in  Table  III.  [Ref.  3] 

Direct  collocation  methods  have  been  used  to  efficiently  solve  low  thrust  orbit 
transfer,  commercial  aircraft  mission  analysis,  chemical  process  control  and  robotics 
problems.  The  Optimal  Trajectories  by  Implicit  Simulation  (OTIS)  program  was  de¬ 
veloped  for  NASA  and  the  U.  S.  Air  Force,  and  has  been  widely  distributed  [Ref.  3]. 
OTIS  uses  a  direct  collocation  method  with  Simpson’s  quadrature  rule  and  Hermite 
cubic  polynomial  approximation  of  the  state  equations  [Ref.  8].  The  OTIS  library  im¬ 
plements  a  sparse  nonlinear  programming  algorithm  using  NPSOL  as  the  NLP  solver. 
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1.  Select  number  of  discretization  points,  k. 

2.  Evaluate  constraints  at  each  node. 

3.  Calculate  the  integral  residual  for  each  step. 

4.  Compute  initial  conditions  ifto- 

5.  Solve  the  resulting  NLP. 

Table  III.  Direct  Collocation  Algorithm 

In  addition,  the  Advanced  Launch  Trajectory  Optimization  Software  (ALTOS)  pro¬ 
gram  developed  for  the  European  Space  Agency  uses  many  of  the  same  features  as 
OTIS. 

C.  DIFFERENTIAL  INCLUSION 

Differential  inclusion  is  very  effective  on  problems  that  have  linear  controls. 
Like  collocation,  differential  inclusion  uses  implicit  integration  rules  to  formulate  non¬ 
linear  constraint  equations  that  are  used  in  the  NLP  formulation  [Ref.  8].  Unlike  col¬ 
location,  however,  differential  inclusion  has  only  been  successfully  implemented  with 
Euler’s  method.  The  linear  controls  allow  the  controls  to  be  represented  explicitly 
in  terms  of  the  states  and  their  rate  of  change  (derivative).  The  method  eliminates 
the  controls  from  the  system  of  equations  by  explicitly  solving  for  the  controls,  and 
then  discretizing  only  the  state  equations  to  reduce  the  number  of  NLP  variables. 
This  is  effective  only  when  an  explicit  formula  for  the  controls  in  terms  of  the  states 
can  be  found.  Since  in  NLP  problems,  the  central  processing  unit  (CPU)  usage  in¬ 
creases  geometrically  with  the  number  of  variables,  eliminating  the  controls  simplifies 
the  problem  and  speeds  convergence  to  a  solution.  This  reduced  problem  size  is  the 
advantage  to  differential  inclusion. 

Low  order  quadrature  rules  are  used  to  integrate  the  state  equations.  Euler’s 
rule,  the  most  common  quadrature  rule  that  has  been  implemented  with  differential 
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inclusion,  has  local  error  0(61 2 3 4 5 6).  Low  order  rules  are  required  to  effectively  isolate 
the  state  and  system  variables  at  each  node  as  a  function  of  only  the  discrete  states. 
The  low  order  of  the  integration  rules  used  in  this  method  require  selection  of  an 
increasingly  smaller  step  size  to  increase  accuracy.  While  decreasing  the  local  error, 
decreasing  the  step  size  also  increases  the  total  number  of  integration  steps  required, 
as  one  step  is  required  for  each  node. 

Another  disadvantage  is  that  if  the  control  histories  are  required,  they  must  be 
solved  for  after  the  optimization.  Differential  inclusion  also  requires  gradient  informa¬ 
tion  during  the  optimization  process  and  requires  prior  knowledge  of  arc  inequalities. 
Prior  knowledge  of  the  arc  inequalities  limits  the  robustness  of  this  approach  and 

calculating  the  gradients  analytically  for  differential  inclusion  can  be  complicated 
[Ref.  8]. 

An  algorithm  using  differential  inclusion  is  shown  in  Table  IV  [Ref.  7],  One 
effective  use  of  differential  inclusion  is  to  attempt  to  solve  a  problem  for  which  collo¬ 
cation  has  failed  to  converge.  The  decreased  problem  size  may  allow  the  new  NLP  to 
converge  [Ref.  8]. 

1.  Solve  for  controls  explicitly  in  terms  of  states  and  their  derivatives. 

2.  Eliminate  control  variables  from  the  system. 

3.  Discretize  the  resulting  constraints  and  inequalities  at  the  nodes. 

4.  Compute  initial  conditions,  ip0. 

5.  Solve  the  resulting  NLP. 

6.  Solve  for  the  resulting  control  histories  (if  required). 

Table  IV.  An  Algorithm  using  Differential  Inclusion 
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D.  SUMMARY  OF  DIRECT  METHODS 

Each  of  the  direct  methods  contain  a  mix  of  advantages  and  disadvantages 
based  upon  their  discretization  and  solution  methodology.  All  direct  methods  share 
the  advantage  of  not  requiring  the  solution  of,  or  solving  for  the  optimality  condi¬ 
tions.  Direct  shooting  and  differential  inclusion  both  tend  to  result  in  smaller  NLPs 
than  those  derived  through  collocation.  Direct  Shooting  is  very  efficient  when  simple 
mathematical  models  suffice  to  control  the  system.  Differential  inclusion  works  very 
well  when  the  controls  are  linear  and  can  be  solved  for  explicitly  in  terms  of  the  state 
equations  and  their  rate  of  change. 

Collocation  effectively  takes  advantage  of  high  order  integration  rules.  The 
advantages  of  high  order  integration  rules  allowing  a  larger  step  size  with  high  accu¬ 
racy  make  collocation  very  robust  and  effective.  The  combination  of  the  accuracy  of 
discretizing  both  the  states  and  controls,  and  the  higher  order  rules  available,  often 
outweigh  the  larger  NLP  formulations,  especially  for  complicated  systems.  We  chose 
to  use  a  collocation  method  in  this  thesis. 
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IV.  A  PSEUDOSPECTRAL  METHOD 


The  Legendre-Gauss-Lobatto  Pseudospectral  method  (LGLP)  is  a  direct  col¬ 
location  method  that  incorporates  spectral  methods  adapted  from  work  in  Computa¬ 
tional  Fluid  Dynamics  (CFD)  [Ref.  9].  LGLP  is  a  direct  collocation  method  that  uses 
a  high  order  integration  rule,  Legendre-Gauss-Lobatto,  and  discretizes  the  time  his¬ 
tories  for  both  the  states  and  controls.  The  use  a  of  high  order  integration  rule  helps 
to  limit  the  number  of  variables  required  for  the  NLP  to  achieve  sufficient  accuracy 
by  allowing  a  large  step  size. 

Pseudospectral  implies  that  the  method  has  spectral  accuracy.  Spectral  ac¬ 
curacy  refers  to  the  property  of  orthogonal  series,  especially  Fourier  Series,  that  the 
kth  coefficient  of  the  expansion  decays  faster  than  any  inverse  power  of  k  when  the 
function  is  smooth  and  the  derivatives  are  periodic.  This  gives  the  Fourier  Series,  and 
other  orthogonal  expansions,  the  property  that  the  truncated  series,  with  a  few  more 
terms,  gives  an  exceedingly  good  approximation  of  the  function,  assuming  proper 
smoothness  of  the  function.  [Ref.  9] 

We  use  orthogonal  Legendre  polynomials  to  approximate  the  state  and  control 
variables.  One  advantage  of  using  orthogonal  polynomials  is  their  close  relationship  to 
Legendre-Gauss-Lobatto  integration  rules.  This  is  exploited  with  the  integral  portion 
of  the  cost  function  and  in  the  implicit  integration  of  the  state  equations  to  transform 
the  optimal  control  problem  into  a  system  of  algebraic  equations.  Polynomial  approx¬ 
imations  of  the  state  and  control  variables  are  used  where  Lagrange  polynomials  are 
the  trial  functions  and  the  coefficients  are  the  values  of  the  state  and  control  variables 
at  the  Legendre-Gauss-Lobatto  (LGL)  points.  The  cost  function  is  discretized  using 
Legendre-Gauss-Lobatto  quadrature  rules. 
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A.  DISCRETIZING  THE  PROBLEM 


The  LGLP  method  can  be  used  for  solving  optimal  control  problems  formu¬ 
lated  as  in  Chapter  II.  The  polynomial  approximations  for  the  state  and  control  func¬ 
tions  are  calculated  in  terms  of  their  values  at  the  LGL  points  lying  in  the  interval 
[-1,1].  We  use  a  transformation  to  express  the  problem  for  r  €  [-1, 1]: 

f  _  ~  *o)r  +  (*/  +  <o)  ,  , 


which  is  equivalent  to 


r  = 


2 1  —  ( tj  +  io) 


(4.2) 


tf  to 

As  a  result  of  the  transformation  t-lrwe  must  adopt  new  symbols  for  the  variables 


x,u,  A,  n,v 

and  the  maps 

Using  the  following  changes  of  variables 

x  y,  u  — »•  v,  A  — »  A,  fi  — >  £,  u  -4  »7, 
and  the  following  change  of  mappings 

J(t)  ->  J(r),£(<)  L(r),  M(t)  ->  7W(r) 

f(t)  ->  <F(r),^0(t)  ->•  SP0(T)iX})f(t)  ^(t). 

It  follows  that  by  using  Equation  (4.1),  Equations  (2.5-2.9)  and  (2.16)can  be  replaced 


II 

AI[y(l),</]+  -  J  L[y(r)Mr)]dr 

(4.3) 

y(r)  = 

('  2t0)lHy(r)Mr))}, 

(4.4) 

!Po(y(-i),to)  = 

0 

(4.5) 

a>/(y(i),*/)  = 

0 

(4.6) 

S(v(t))  < 

0 

(4.7) 

ft(x,  A,v)  = 

(tzjf£)yl^+(f/-to')L  +  <Tg 

(4.8) 
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B.  LEGENDRE  POLYNOMIALS 

There  are  three  main  classes  of  classic  orthogonal  polynomials,  Jacobi,  Laguere 
and  Hermite  polynomials  [Ref.  10].  Legendre  and  Chebyshev  polynomials,  both  used 
in  collocation  schemes  for  direct  methods,  are  of  the  Jacobi  polynomials  class.  One 
advantage  of  using  orthogonal  polynomials  is  that  their  derivatives  also  form  a  set 
of  orthogonal  polynomials  [Ref.  10].  The  polynomials  are  bounded  in  every  proper 
subinterval  of  (—1, 1),  which  allows  a  transformation  to  Z-2[ — 1, 1],  to  ensure  the  states 
and  controls  will  be  bounded. 

The  Rodriguez  formula  is  commonly  used  to  express  Jacobi  polynomials 

(1  - 1)0(1  +  xypv^z)  =  [(!  -  +  *)n+1r]  10]-  (4.9) 

Pn(x )  represents  the  orthogonal  polynomial  of  the  nth  degree  evaluated  at  the  point 
x.  The  variables  (3  and  7  are  parameters  which  define  the  type  of  polynomial  defined 
by  the  Rodriguez  formula.  When  /?,  7  >  —1,  the  polynomials  expressed  by  the  for¬ 
mula  are  Jacobi  polynomials.  The  weight  function  that  characterizes  the  orthogonal 
polynomial,  (1  —  a:)/3(  1  +  a:)7,  is  also  incorporated  in  Equation  (4.9).  The  Legen¬ 
dre  Polynomials  are  defined  by  j3  =  7  =  0  resulting  in  a  weight  function  w(x)  =  1 
[Ref.  11].  The  Rodriguez  formula  for  Legendre  polynomials  simplifies  to 

(4.10) 

Thus,  Ln(x)  is  the  Legendre  polynomial  of  degree  n  with  leading  term  The 

Legendre  Polynomials  may  be  generated  by  the  three  term  recurrence  relationship 

Ob  4. 1  u 

Lk+i(x)  =  -j— j ~xLk(x)  -  j-^-Lk-i(x)  (4.11) 

with 

Lo(x)  =  1  and  L\{x)  —  x. 

The  orthogonal  polynomials  have  a  close  relationship  with  the  theory  of  real 
tridiagonal  matrices.  This  relationship  allows  new  numerical  linear  algebra  codes  that 


25 


exploit  the  matrix  structure  to  find  the  zeros  of  the  Legendre  polynomials.  The  roots 
Xi  ?  i  =  1, 2,  —  ,  n  are  the  eigenvalues  of  the  tridiagonal  matrix 


J n  — 


$1  72 

72  s2  •  • . 

In 

7n 


(4.12) 


The  unique  orthogonal  properties  of  the  Legendre  polynomials  allow  us  to 
exploit  their  boundedness  and  structure  by  transforming  the  original  problem  space 
into  [—1,1].  The  Legendre  polynomials  satisfy  the  ODE 


(1  -  x2)Ln  -  2 xLn  +  n(n  +  l)y  =  0  [Ref.  11]. 
With  Legendre  polynomials  we  have  the  properties  that 

I-i  x)Ln(x)dx  =  0  m  n 

/x  2 

Lm(x)Ln(x)dx  =  - - - 

-i  zn  +  1 


m  —  n. 


(4.13) 


(4.14) 


(4.15) 


This  allows  us  to  neglect  cross  terms  in  the  expansions  of  the  state  and  control  vari¬ 
ables.  When  a  function  is  expanded  with  Legendre  polynomials,  the  approximation 
formulas  are 

OO 

/(*)  =  H  anLn(x)  (4.16) 

n=0 

and 

2n  +  1  /i 


dn  — 


■J  f(x)Ln(x)dx. 


(4.17) 


C.  APPROXIMATING  THE  STATES  AND  CONTROLS 

Let  Ln(t)  be  the  Legendre  polynomial  of  degree  N  on  the  interval  [-1,1]. 
In  the  Legendre  collocation  approximation  of  (4.3)-(4.8),  we  use  the  LGL  points, 
Th  l  =  0, 1, . . . ,  N  which  are  given  by 


T0  =  -1 


26 


Ti,  r2,  .  .  .  Tjv-l, 


the  zeros  of  Ln,  the  derivative  of  the  Legendre  polynomial  Ln,  and 

tn  =  1. 

For  approximating  the  continuous  equations,  we  seek  polynomial  approximations  for 
the  state  and  control  equations.  We  define 


yN(T)  =  £y 

l- 0 

v"(r)  =  y  v(r«)*!TI, 


(4.18) 


(4.19) 


where,  for  /  =  0, 1, ...,  N, 

“■udiss^  <“> 

are  the  Lagrange  polynomials  of  order  N.  Notice  that  the  Lagrange  polynomials 


(4.20) 


require  our  previous  calculation  of  Ln  and  Ln.  It  can  be  shown  that 


4>l{Tk)  =  &lk  = 


1  if  l  =  k 


(0  if  /  ^  k. 

From  this  orthogonal  property  of  (j>i  it  follows  that 

y  N(n)  =  y(ri) 

vn(ti)  =  v(r<). 


(4.21) 

(4.22) 


Thus  the  value  of  the  interpolating  polynomial  at  the  LGL  points  is  exactly  the  value 
of  the  continuous  state  and  control  functions  evaluated  at  the  LGL  points. 

To  facilitate  the  NLP  formulation,  we  use  the  notation 

a;  :=  y(r,),  b,  :=  v(r;), 


to  rewrite  (4.18)-(4.19)  in  the  form: 


yN(r)  = 

i—o 

vJV(r)  =  J2 


(4.23) 


(4.24) 
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D.  CALCULATING  DERIVATIVES 

To  express  the  derivative  y n(t)  in  terms  of  y n(t)  at  the  collocation  points 
T/,  we  differentiate  (4.18).  This  differentiation  of  the  approximating  polynomial  is 
accomplished  by  a  matrix  multiplication  of  the  following  form: 

yN(n)  =  32  Dkiy(n),  (4.25) 

/=0 

where  Dki  are  entries  of  the  (N  +  1)  x  (N  +  1)  differentiation  matrix  D 


An  example  of  the  Legendre  polynomial  and  the  differentiation  matrix  for  a 
fixed  N  is  illustrative  of  this  technique.  Letting  N  =  4,  we  have  the  resulting  4th 
order  Legendre  polynomial, 


L4  =  -(35a:4  -  30a:2  +  3). 

O 


(4.27) 


Dividing  the  interval  into  four  time  periods  and  solving  for  the  Lobatto  points,  the 
zeros  of  L4 ,  results  in  r0  =  -1,  tx  =  -.6547,  r2  =  0,  r3  =  0.6547,  r4  =  1.  The 
differentiation  matrix  results  in  a  matrix  of  the  following  form 

-5  6.7565  -2.6667  1.4102  -0.5 

-1-2410  0  1.7457  -0.7638  .2590 

0.3750  -1.3366  0  1.3366  -0.3750  .  (4.28) 

-0.2590  0.7638  -1.7457  0  1.2410 

0.5  -1.4102  2.6667  -6.7565  5 


28 


E.  DISCRETIZING  INTEGRALS 

Next,  the  integral  (4.3)  is  discretized.  Substituting  (4.23)  and  (4.24)  in  (4.3) 
and  using  the  Gauss-Lobatto  integration  rule,  we  obtain 

JN(a,b,t/)  =  M(yN(l),tf)  +  L(y N,vN)dr 

=  M{yN(l),tf)  +  tf  -  -  J2  L  (lla iMn),  wk 

=  M(aN,tf)  +  —  -  Y,M*k,bk)wk  (4.30) 

1  k= o 

The  last  equality  is  obtained  from  4>i(tk)  =  $ik-  The  coefficients  are  a  =  (ao,  ai, . . . ,  a^),and 
b  =  (bo,  bj, . . . ,  b^).  The  weights  are  given  by 


N(N  + 1)  [Lis[{Tk))2 


0,1,...,  N. 


F.  THE  PARAMETRIC  OPTIMIZATION  PROBLEM 

The  state  equations  and  the  initial  and  terminal  state  conditions  are  discretized 
by  first  substituting  (4.23)-(4.24)  in  (4.8)  and  collocating  at  the  LGL  points,  r*.  Using 
the  notation  for  a  and  b,  the  state  equations  are  transformed  into  the  following 
algebraic  equations 

Afc(a,b)  =  ^y-^F(afc,bfc)  -  ck  =  0,  k  =  0,1,..., iV,  (4.31) 
where  ck  is  as  defined  in  (4.29),  and  the  initial  conditions  are 

®ro(yAr(-l),<o)  =  0  or  (4.32) 

?Po(ao,*o)  =  0.  (4.33) 
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The  terminal  state  conditions  are 


S>/(y"(l),</)  =  0  or  (4.34) 

tj)  =  0  (435) 

The  control  inequality  constraints  are  approximated  by 

g(vN(r,))  <  0,  k  =  0,l,...,N,  or  (4.36) 

g(b*)  <  0,  k  =  0,1,..., N.  (4.37) 


To  summarize,  the  optimal  control  problem  (4.3)-(4.8)  is  approximated  by  the 


following  nonlinear  optimization  problem:  Find  the  coefficients 

a=  (ao,ai,...,ajv)  (4.38) 

b  =  (b0,  bi, . . . ,  b/v)  (4.39) 

and  possibly  the  final  time  tj ,  to  minimize  the  cost  function 

J"(a,  b)  =  A4(aw ,  t,)  +  £  L(atj  hk)vik  (4 .40) 

Z  k=0 

subject  to 

Afc(a,  b)  =  (iL_io)jF(afc,bfc)-cA  =  0,  k  =  0,1,...,  JV,  (4.41) 

B*(b)  =  g(bfc)  <  0,  k  =  0,1,..., N,  (4.42) 

&o(ao,to)  =  0,  (4.43) 

9f(a.N,tf)  =  0.  (4.44) 


The  LGLP  method  relies  upon  a  simple  conversion  that  maintains  much  of 
the  structure  of  the  original  problem.  By  collocating  at  the  LGL  points  the  functions 
are  evaluated  without  any  dependence  on  neighboring  points.  Once  the  conversion 
to  a  parameter  optimization  problem  is  complete,  the  system  may  be  solved  by  using 
an  NLP  solver. 
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V 


THE  LGLP  GRAPHICAL  USER 
INTERFACE 

The  implementation  of  the  LGLP  algorithm  in  MATLAB,  using  the  opti¬ 
mization  toolbox,  requires  a  function  file  containing  the  state  equations,  constraints 
and  cost  function  as  input.  Each  state  equation  must  be  transformed  and  recorded 
in  the  input  file.  This  requires  a  detailed  knowledge  of  the  algorithm  to  compute 
the  transformation  to  [—1,1],  a  familiarity  with  MATLAB  function  construction,  file 
input /output  (I/O),  and  constr.mfrom  the  optimization  toolbox.  Although  a  MAT¬ 
LAB  knowledge  base  may  be  assumed,  and  the  published  algorithm  may  be  learned, 
we  desire  to  provide  an  interface  that  allows  a  user  to  input  information  about  an 
optimal  control  problem  in  the  standard  form  and  have  the  program  transform  and 
solve  the  problem. 

A.  MATLAB  GUIS 

GUI  design  has  almost  become  synonymous  with  windows  programming.  With 
programming  for  Windows,  UNIX  or  Apple,  the  operating  system  provides  the  user 
with  a  familiar  point  and  click  interface  to  most  system  functions.  Most  operating 
systems  (OS)  provide  access  to  a  command  line  interface,  but  as  each  succeeding  gen¬ 
eration  of  computers  increasingly  relies  upon  GUIs,  portability  to  multiple  platforms 
and  accessibility  to  wide  audiences  will  only  be  achieved  through  GUIs. 

MATLAB  at  its  core  provides  a  command  line  interface  to  allow  a  C-like 
code  to  be  easily  written  and  quickly  interpreted  to  solve  mathematical  problems, 
especially  in  the  fields  requiring  numerical  linear  algebra.  The  creation  of  a  GUI 
within  MATLAB  allows  access  to  the  strength  of  the  command  line  interface  while 
providing  a  dynamic  interface  to  the  solution  of  a  problem. 

MATLAB  GUI  design,  and  GUI  design  in  general,  center  on  several  principles 
presented  in  the  MATLAB  Graphics  and  GUIs  manual  [Ref.  15].  A  GUI  should  have 
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the  characteristics  shown  in  Table  V. 


1.  Simple 

2.  Consistent 

3.  Familiar 

4.  Dynamic 


Table  V.  GUI  Design  Characteristics 

The  GUI  must  be  simple  to  use,  as  the  goal  is  to  create  an  interface  within 
MATLAB  which  makes  solving  a  problem  easier  than  using  the  command  line  inter¬ 
face.  The  project  must  be  consistent  to  prevent  confusion  as  the  problem  is  solved. 
A  familiar  interface  used  throughout  the  GUI  will  accelerate  the  learning  curve,  al¬ 
lowing  a  user  to  quickly  accomplish  new  tasks  by  mimicking  completed  ones.  Lastly, 
the  interface  must  be  dynamic,  to  allow  the  user  to  efficiently  maneuver  through  the 
solution  process,  changing  incorrect  values,  and  providing  feedback  for  actions  while 
solving  the  problem.  A  GUI  must  also  allow  the  user  to  directly  access  the  input  files. 
Direct  file  access  is  required  for  advanced  users  who  might  not  want  to  wade  through 
several  layers  of  point  and  click  interface  to  change  a  single  variable  or  parameter 
before  rerunning  an  optimization. 

The  hierarchy  of  control  for  MATLAB  graphics  provides  the  starting  point  for 
manipulation  of  GUIs,  Figure  1.  Each  object,  figure,  user  interface  control  (uicontrol), 
axes,  and  user  interface  menu  (uimenu)  have  properties  and  callbacks  that  control 
its  appearance  and  function.  The  GUIDE,  Figure  2,  simplifies  access  to  the  handle 
graphics  objects  and  allows  creation  of  a  GUI  by  clicking  and  dragging  the  desired 
components  off  of  a  visual  menu.  The  GUIDE  control  panel  provides  the  building 
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Screen/Root 

Uicontrols 

Axes 

Uimenus 

Figure  1.  Hierarchy  of  MATLAB  Graphics 

blocks  of  the  GUI  and  allows  control  of  the  current  figures.  Figures  are  the  objects 
that  contain  all  other  MATLAB  graphics  objects. 

Each  of  the  objects  will  have  a  unique  graphics  handle,  represented  by  an 
integer  or  a  float.  Before  an  object  is  manipulated,  it  must  be  grabbed,  much  like 
a  cup  must  be  grabbed  before  its  contents  may  be  drunk.  Several  reserved  words  in 
MATLAB  make  manipulating  the  handle  graphics  very  simple.  The  abbreviations 
gcf,  gca,  and  geo  are  used  to  refer  to  the  active  figure,  axes  and  object  respectively. 
The  object  that  initiates  a  callback  is  referred  to  as  gebo.  In  windows  programming, 
handles  are  used  to  control  objects.  To  agree  with  convention,  graphics  handles  are 
referred  to  as  handle  graphics.  Handle  graphics  are  covered  in  detail  in  the  MATLAB 
documentation  and  In  Graphics  and  GUIs  with  MATLAB  [Ref.  13]. 

The  alignment  tool  allows  easy  alignment  of  the  controls  of  the  GUI.  The 
menu  editor  provides  access  to  the  figure’s  user  defined  menus,  and  the  callback 
editor  controls  the  functions  that  will  be  executed  by  each  of  the  figures  and  the 
GUI  controls.  The  property  editor  provides  access  to  every  property  of  the  figure, 
its  uicontrols  and  its  axes.  The  property  editor  accesses  each  of  the  properties  of  the 
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Figure  2.  MATLAB  Guide 


objects.  Their  size,  color,  and  functionality  can  all  be  edited  from  within  the  property 
editor. 

The  callback  editor  is  the  tool  that  creates  the  true  functionality  for  the  GUI. 
The  callbacks  associated  with  each  object  are  shown  in.  Table  VI. 

The  create  function  works  efficiently  to  provide  default  values  for  variables, 
editable  fields,  and  other  uicontrols.  Each  time  a  figure  or  handle  graphic  is  created, 
the  create  callback  function  is  executed.  A  simple  create  function  call  from  a  text 
edit  box  of  set  Cgcbo,  ’String1  ,Names_States)  allows  a  variable,  Names_States,  to 
be  written  into  the  text  edit  box  that  had  just  been  created  and  called  the  create 
callback. 
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1.  Callback  Function 

2.  Buttondown  Function 

3.  Create  Function 

4.  Delete  Function 


Table  VI.  GUI  Callbacks 

This  usage  of  the  create  function  places  the  current  states  names  in  the  window 
so  the  user  can  provide  initial  guesses  for  the  states  and  control  variables.  This 
enhances  the  dynamic  nature  of  the  interface  by  incorporating  current  values  and  not 
requiring  redundant  input  of  values  while  still  clearly  illustrating  the  current  values. 

The  callback  function  executes  when  a  button  is  pushed,  a  slider  slid,  a  box 
checked  or  field  edited.  A  simple  callback  is  illustrated  by  the  State  Equations  button 
on  the  Function  File  GUI,  Figure  5.  The  code  StatesEquationsFile  is  executed 
and  the  State  Equations  figure  is  displayed.  Another  call  back ,  close(gcbo),  closes 
the  State  Equations  box  and  returns  to  the  Function  File  GUI. 

The  delete  function  may  be  used  to  clear  the  workspace  when  an  object  is 
closed.  This  can  keep  the  workspace  from  becoming  cluttered  with  variables  that  are 
no  longer  needed.  The  buttondown  function  activates  code  when  the  cursor  is  located 
over  the  edge  of  any  object.  This  differs  from  the  callback  function,  as  the  mouse 
does  not  click  on  the  object,  but  only  needs  to  rest  over  the  object  to  activate  the 
function  call.  These  basic  tools  of  the  GUIDE  allow  for  the  rapid  development  of 
a  GUI  and  provide  a  framework  that  can  be  edited  by  hand  to  further  optimize  or 
tailor  the  code. 

Several  usability  features  and  design  considerations  are  necessary  when  cre¬ 
ating  large  MATLAB  GUIs.  Createfunctions,  callback  functions,  and  file  operations 
form  the  backbone  of  the  GUI.  Designing  the  callbacks  involves  deciding  to  use  a 
script,  function,  or  in-line  code  for  all  the  executable  code.  In-line  code  has  the 
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disadvantage  that  long  sections  of  code  are  very  difficult  to  read  in  the  GUIDE  en¬ 
vironment  and  are  not  stored  in  the  GUI  m-file  (ASCII),  but  instead  in  the  mat-file 
(binary).  This  makes  the  code  editable  only  within  the  GUIDE1  environment  and  pre¬ 
vents  editing  the  code  by  hand.  Function  calls  execute  code  the  fastest,  but  require 
proper  parameter  passing  to  ensure  the  proper  variables  are  passed  into  the  function 
workspace.  A  function  space  is  created  each  time  a  function  is  executed.  This  pro¬ 
vides  safety  from  overwriting  variables  in  the  main  workspace,  and  ensure  that  all  the 
necessary  variables  will  be  available  when  the  function  is  executed.  Writing  individual 
functions  can  require  a  large  number  of  files  to  be  associated  with  the  GUI.  Scripts 
are  slower  than  functions,  but  by  using  a  switch  statement  in  a  script,  numerous 
callbacks  can  be  made  through  a  single  file  to  reduce  the  number  of  files  included  in 
the  implementation.  A  switch  statement  could  also  be  used  with  a  function  call,  but 
the  parameter  passing  becomes  increasingly  difficult  as  different  calls  are  made  to 
a  function  with  different  parameters.  MATLAB  does  not  have  strict  type  checking, 
so  generic  parameters  may  be  used  with  a  function,  but  the  code  is  very  cryptic  as 
the  variable  names  will  not  be  named  according  to  their  function.  For  example,  a 
MATLAB  function  could  take  three  parameters, 

function  switcherCparal,  para2,  type) ; 

and  have  a  switch  statement  within  the  function  that  evaluated  the  value  of  the 
parameter  type.  For  one  value  of  type,  the  other  two  parameters,  paral ,  para2, 
could  be  strings.  For  another  value,  the  parameters  could  be  floats.  Neither  of  the 
parameter  names,  paral ,  para2,  provide  any  insight  into  their  values. 

B.  GUI  DESIGN 

The  creation  of  a  GUI  for  the  LGLP  algorithm  makes  the  method  accessible 
to  a  wider  audience  of  engineers  and  applied  mathematicians.  The  mathematics 
required  to  perform  the  transformation  of  the  problem  into  [-1,1],  the  derivation 
and  application  of  the  differentiation  matrix,  and  the  implicit  integration  all  provide 
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potential  stumbling  blocks  to  easy  access  to  the  strength  of  the  LGLP  algorithm. 
These  stumbling  blocks  are  removed  by  an  efficient  GUI. 

A  GUI  must  have  unity  of  purpose.  The  purpose  of  this  interface  is  to  solve 
optimal  control  problems.  The  interface  must  be  simple  and  consistent  to  allow  easy 
solution  of  problems  and  prevent  the  interface  from  distracting  from  the  solution 
process.  The  original  implementation  required  writing  a  function  file  outlining  the 
states,  constraints  and  cost  functions,  and  then  creating  a  function  or  script  to  set 
the  proper  variables  in  the  main  workspace  and  invoke  the  NLP  solver.  For  simplicity 
we  kept  this  two  step  process  in  the  design  of  the  GUI,  Figure  3. 

The  first  design  task  for  the  new  GUI  for  the  LGLP  algorithm  was  to  sketch 
the  interface.  The  two  main  tasks,  creating  the  file,  and  setting  the  variables,  should 
be  separated  for  clarity.  Also,  once  the  function  file  has  been  created  the  optimization 
window  must  allow  for  repeated  runs  in  multiple  sessions  without  recreating  the  file 
code. 


Figure  3.  Initial  GUI  Design 

Figure  4  serves  as  the  starting  point  for  the  GUI.  Typing  opt  at  the  command 
line  will  execute  the  GUI.  In  Figure  5,  the  Function  File  figure  creates  the  input  files 
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Figure  4.  LGLP  GUI  Welcome  Figure 

necessary  for  NLP  solver,  constr.m  [Ref.  16].  The  GUI  must  be  able  to  convert 
an  optimal  control  problem,  with  fixed  or  free  final  time,  all  of  its  state  equations, 
control  constraints,  equality  and  inequality  constraints,  and  parameters  into  an  NLP. 

Figure  6  solves  the  optimization  problem.  The  optimization  of  the  input  file 
requires  the  parameters  for  constr .  m  to  be  selected.  Also,  the  interval  size  for  the 
discretization  must  be  selected,  the  values  for  any  parameters  must  be  set,  and  the 
initial  guesses  must  be  provided  for  each  of  the  states,  controls,  and  final  time. 

Another  design  decision  was  to  ensure  that  default  values  would  be  created 
to  prevent  errors  if  not  all  of  the  possible  inputs  were  selected.  By  limiting  the 
callbacks  executed  by  the  main  figure’s  children  to  merely  receiving  input,  the  main 
computation  could  be  isolated  to  the  main  GUI  figures.  This  ensured  that  if  all  of 
the  variables  had  been  input,  then  the  GUI  would  run  properly. 

The  main  input  for  the  constr  .m  optimization  algorithm  is  a  file  that  contains 
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Figure  5.  Function  File  Creator  Figure 

the  state  equations,  constraints,  and  cost  function.  This  requires  MATLAB  to  create 
and  manipulate  a  file  through  use  of  file  I/O  that  mirrors  ANSI  C  or  C++.  A  second 
MATLAB  file  feature  that  simplified  the  GUI  was  the  mat-file  constructs  which  saves 
the  values  of  the  variables  used  during  the  optimization.  Every  time  an  m-file  was 
created,  a  sister  mat-file  was  created  with  the  variables  declarations.  This  allowed 
all  of  the  variables  associated  with  a  function  file  to  be  loaded  into  the  workspace 
with  create  functions.  Once  a  filename  has  been  selected,  the  GUI  attempts  to  load 
a  corresponding  mat  file.  A  flag  returned  by  MATLAB  when  the  mat  function  does 
not  exist,  which  occurs  the  first  time  a  function  is  created,  prompts  the  initialization 
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Figure  6.  Optimization  Figure 

of  all  the  variables  required  by  the  GUI.  The  diary  function  was  also  used  to  create 
an  output  file  after  the  optimization  without  requiring  redundant  file  I/O.  The  diary 
file  captures  all  the  optimization  results  that  display  in  the  MATLAB  workspace, 
creating  a  record  of  each  optimization  run.  By  choosing  the  append  option  for  file 
I/O,  each  optimization  run  is  appended  to  the  same  diary  file.  The  clock  function  is 
used  to  record  the  date  and  time  of  day  of  optimization  in  the  diary  file  to  distinguish 
between  runs  in  the  file,  see  Appendix  A,B,  and  C. 

Every  editable  textbox  in  a  GUI  takes  a  string  as  input.  MATLAB  does  not 
allow  a  particular  textbox  to  be  declared  as  a  float.  This  convention  of  all  inputs 
being  declared  as  strings  motivates  a  style  of  coding  that  takes  advantage  of  this 
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structure.  Utility  functions  designed  to  work  with  strings  and  accomplish  simple, 
repetitive  tasks  can  streamline  the  coding  of  a  GUI.  A  multiline  input  box  returns 
a  matrix  of  string  variables.  In  a  GUI  with  two  multiline  textboxes,  we  assign  the 
strings  in  one  box  as  variable  names  and  use  the  inputs  from  the  other  box  for  values 
of  these  variables.  The  function  declare. m  is  a  simple  utility  function  that  takes  two 
matrices  of  strings  and  combines  the  strings  into  code  that  will  assign  the  names  of 
the  variables  to  the  values,  see  Appendix  A.  The  eval  command  is  used  in  MATLAB 
to  execute  the  commands  contained  in  a  string  variable.  The  eval  is  used  to  execute 
the  string  matrices  that  are  returned  from  declare. m. 

Two  other  utility  functions  were  created  with  this  GUI,  optcall  .mand  costcall  .m. 
Both  accept  parameters  as  input  to  create  function  calls  to  be  evaluated  by  the  GUI. 

The  function  opt  call,  m  creates  a  function  call  for  constr.m,  and  costcall  .m  eval¬ 
uates  the  cost  function  for  the  files  created  with  the  File  Creator  GUI. 

Creating  help  files  for  a  GUI  is  an  important  part  of  the  design.  MATLAB 
has  two  features  for  creating  help,  comments  within  an  m-file  and  the  Contents. m 
file.  The  declare. m  code  in  Appendix  A  illustrates  the  MATLAB  comment  help 
convention.  The  initial  comments  in  an  m-file  after  the  function  header  are  included 
by  MATLAB  as  help  for  the  function.  When  help  declare  is  typed  at  the  command 
line  the  initial  commented  lines  of  code  in  the  file  is  displayed.  MATLAB  displays 
the  comments  until  the  first  line  after  the  function  header  that  does  not  include  the 
comment  character. 

The  second  convention  is  the  Contents  .mfile.  The  comments  in  the  Contents  .m 
file,  are  displayed  when  the  help  command  is  called  with  a  directory  as  its  input  argu¬ 
ment.  When  a  directory  does  not  include  a  Contents  .mfile,  the  first  help  line  (Hi)  of 
each  m-file  in  the  directory  is  displayed.  Creating  a  Contents  .m  file  provides  another 
method  for  the  user  to  get  help  on  a  GUI,  see  Appendix  B.  All  GUI  files  may  have 
comments  and  help  lines  built  into  their  code.  However,  when  a  GUI  is  edited  with 
the  GUIDE ,  any  comments  are  erased  and  a  standard  comment  explaining  that  this 
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code  is  a  computer  generated  object  is  inserted.  This  severely  limits  the  usefulness 
of  adding  any  comments  lines  within  the  GUI.  By  placing  the  comments  that  would 
be  included  within  the  GUI  in  the  Contents. m  file,  the  comments  will  be  preserved 
regardless  of  the  number  of  times  a  GUI  is  modified. 

A  simple  example  problem  will  illustrate  the  use  of  the  GUI.  The  GUI  was 
laid  out  to  encourage  all  text  boxes  to  have  values  entered  before  the  various  buttons 
were  pushed.  All  the  buttons  are  designed  to  be  executed  from  top  to  bottom,  with 
a  familiar  Finished  button  in  the  bottom  of  each  figure. 
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VI. 


A  GUI  TUTORIAL-THE  CART 
PROBLEM 

In  a  recent  paper  by  Conway  and  Larson,  “Collocation  Versus  Differential  In¬ 
clusion  in  Direct  Optimization”,  the  following  problem  was  posed  [Ref.  8].  A  cart  is 
placed  on  a  track  initially  at  rest.  An  external  force  u(t)  is  applied  to  the  cart  of  unit 
mass.  The  cart  is  subject  to  drag  depending  linearly  on  the  velocity.  The  system 
equations  are 


ii  =  x2  (6.1) 

x2  =  —x2  +  u.  (6.2) 

The  force  u(f)  is  to  be  applied  so  as  to  satisfy  a  terminal  constraint  at  the  final  time 
t j  which  is  a  combination  of  position  and  velocity 

4>j  =  axi(tj)  +  bx2(tf)  —  c  =  0 

while  minimizing  the  integral  of  the  square  of  the  control 

J  =  f  u2dt.  (6-3) 

Jto 

With  the  values  of  the  constants  chosen  as  a  =  1.0,  b  =  —2.694528,  and  c  —  —1.155356, 
and  tf  =  2,  the  value  of  the  objective  function  is  J  =  0.577678  [Ref.  8].  This  problem 
provides  an  easy  tutorial  as  the  analytic  solution  is  known  and  the  problem  is  quickly 
optimized  by  constr.m. 

A.  FILE  CREATION 

Typing  opt  at  the  command  line  starts  the  GUI  and  allows  access  to  Create 
Optimization  File  button.  Pressing  this  button  starts  the  input  file  creation  process. 
The  implementation  begins  with  the  creation  of  a  file  that  contains  the  constraints, 
state  equations,  and  cost  function  for  the  problem.  The  first  step  is  always  to  enter  a 
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filename  for  the  file;  see  Figure  7.  Entering  the  filename  will  set  the  current  filename 
and  will  also  check  to  see  if  a  MATLAB  data  file  exists  for  this  filename.  If  a  file 
exists,  the  file  will  be  loaded.  This  loads  all  the  variables  associated  with  the  problem 
into  the  MATLAB  workspace. 


Figure  7.  Cart  Problem  File  Creator  -  GUI 


After  the  filename  for  the  input  file  to  constr.m  has  been  chosen,  buttons 
execute  additional  figures  to  collect  the  necessary  variables  and  parameters  for  the 
problem.  Next  the  state  variables,  and  x2,  are  entered  into  the  GUI  implementation 
of  our  method.  Pressing  the  State  Equations  button  will  open  a  window  for  inputting 
the  state  equations  for  aq  and  X2,  see  Figure  8.  The  state  equations  are  entered  with 
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Figure  8.  Cart  Problem  State  Equations  -  GUI 

the  x  assumed  as  the  left  hand  side  of  the  equality.  Here  we  must  make  a  transition 
into  MATLAB  numerical  notation.  All  entries  must  be  legal  MATLAB  expressions. 
Of  special  note  is  the  difference  between  matrix  multiplication  and  element-by-element 
array  multiplication  or  dot  multiplication.  Dot  multiplication  requires  the  use  of  the 
.*  construct  to  ensure  dimensional  agreement  within  MATLAB.  The  states,  X\  and  x2, 
are  both  n  x  1  vectors.  Each  component  corresponds  to  a  discrete  time  step  between 
the  initial  and  final  times.  A  simple  system  of  equations  arises  from  choosing  to  use 
three  LGL  points  for  this  example.  The  state  equations,  (6.1  -6.2),  can  be  expressed 
as  a  systems  of  six  equations  of  scalar  variables: 


Zi(l) 

=  *2(1) 

(6.4) 

ii(2) 

=  *2(2) 

(6.5) 

ii(3) 

=  *2(3) 

(6.6) 

*2(1) 

=  -*2(1)+«(1) 

(6.7) 
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•£2(2)  —  —£2(2)  +  w(2) 
*2(3)  =  -*2(3)  +  ti(3). 


(6.8) 

(6.9) 


This  form  of  the  systems  of  equations  illustrates  that  the  operations  are  con¬ 
ducted  component-wise.  A  slightly  different  equation  illustrates  dot  multiplication. 


Let  the  second  state  equation  be 

x2  =  -x2xi  +  u.  (6.10) 

The  second  state  equation  can  be  expressed  as 

*2(1)  =  -*2(l)sj(l)  +  u(l)  (6.11) 

*2(2)  =  — ®2(2)a:1(2)  +  u(2)  (6.12) 

*2(3)  =  -ar2(3)®1(3)  +  u(3).  (6.13) 


The  operations  are  component-wise  in  this  system  of  equations.  To  properly  represent 
this  system  in  MATLAB  we  use  dot  multiplication,.*.  The  proper  MATLAB  input 
for  the  GUI  would  be  -x2(n)  .*xl(n)u(n)  .  MATLAB’s  dot  multiplication  performs 
the  necessary  component-wise  operations  on  the  state  and  control  vectors  without 
requiring  that  the  continuous  equations  be  converted  into  a  system  of  discrete  equa- 
tions. 

When  using  the  Windows  OS,  a  CTRL  ENTER  must  be  used  to  record  a 
multi  line  input.  MATLAB  only  distinguishes  between  single  and  multi-line  inputs 
and  the  CTRL  ENTER  tells  MATLAB  the  user  has  finished  with  the  text  edit  box. 
Extra  lines  should  not  be  added  between  the  values  or  after  the  values.  The  UNIX 
OS  does  not  require  a  CTRL  ENTER,  but  if  a  standard  ENTER  is  used  after  the 
last  entry,  and  extra  empty  entry  will  be  added  to  the  input  vector. 

After  closing  the  state  equations  window,  the  Controls  variable  button  is 
pressed  and  the  control  u  is  entered  in  the  Control  Variables  window,  see  Figure 
9.  To  set  the  variables  x\  and  x2  to  the  initial  value  of  xi  =  0  and  x2  —  0,  the  first 
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Figure  9.  Cart  Problem  Controls  -  GUI 

discrete  value  of  the  X\  and  x2  vectors  must  be  set  to  zero.  This  is  accomplished  by 
writing  the  equation  in  a  standard  NLP  form  with  the  equation  set  equal  to  zero. 

xi  —  0  =  xl(l)  (6.14) 

x2  =  0  =  x2(l).  (6.15) 

For  a  slightly  different  initial  condition,  x\  =  1,  we  enter 

xl(l)  -  1.  (6.16) 

It  is  assumed  (6.16)  is  set  equal  to  zero.  This  format  is  standard  for  NLP  formulations 

and  will  be  used  for  all  the  constraints  entered  into  the  GUI.  The  two  initial  conditions 
for  the  cart  problem,  Xi  =  0  and  x2  —  0,  are  entered  in  the  Initial  Conditions  window 
in  Figure  10. 

Next,  the  final  time  boundary  conditions  are  entered  along  with  a  selection 
of  fixed  or  free  final  time,  see  Figure  11.  The  discretized  final  time  variables  are 
represented  as  the  last  component  of  the  state  vectors,  in  this  example  xl(n)  and 
x2(n).  Within  the  GUI  final  time  constraints  and  initial  time  constraints  are  both 
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Figure  10.  Cart  Problem  Initial  Conditions  -  GUI 

treated  as  scalar  constraints.  The  boundary  conditions  are  input  in  two  figures  for 
clarity  only.  The  difference  between  boundary  conditions  is  transparent  to  the  NLP. 

Next,  the  three  parameters  a,  6,  and  c  are  entered  in  the  parameters  window, 
see  Figure  12.  Values  needed  by  the  function  while  is  it  being  optimized  must  be 
passed  through  the  parameters  window.  The  values  of  the  parameters  are  free  at  the 
time  of  the  file  creation.  The  values  are  set  in  the  Optimization  figure,  allowing  for 
easy  manipulation  of  the  constants  while  solving  the  problem. 

Finally,  the  cost  function  is  entered.  The  type  of  cost  function  is  selected, 
in  this  case  Lagrange.  The  cost  function  requires  raising  the  control  to  the  second 
power,  see  Figure  13.  This  is  accomplished  by  component  wise  multiplication  or 
exponentiation.  As  explained  earlier,  the  equation  requires  the  .  ~  construct  to 
square  the  control.  These  operations  are  done  element-by-element  and  are  expressed 
in  MATLAB  as  u.*2or  u.*u. 

The  selection  of  a  type  of  cost  function  controls  the  number  of  inputs  allowed 
to  the  Cost  Function .  If  either  the  Mayer  or  Lagrange  type  are  entered,  only  one 
input  is  allowed.  When  the  Bolza  type  cost  function  is  selected,  two  inputs  are 
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Figure  11.  Cart  Problem  Final  Conditions  -  GUI 


Figure  12.  Cart  Problem  Parameters  -  GUI 


expected.  With  a  Bolza  cost  function,  the  first  line  input  is  the  integral  portion  of 
the  cost  function,  and  the  second  line  of  the  input  is  the  scalar  portion.  The  proper 
operation,  addition  or  subtraction,  must  be  included  before  the  scalar  portion  of  a 

Bolza  cost  function.  The  proper  input  of  the  Bolza  cost  function  ///  u2  dt  +  1  would 
be 

u.~2 

+  1. 

For  the  Mayer  cost  function  tj  the  equivalent  Lagrange  cost  function  would  be 
ones(l,n)  to  allow  for  the  numeric  integration  from  0  to  tf.  The  MATLAB  code 
ones(l  ,n)  creates  an  n  dimensional  vector  with  each  component  equal  to  one.  The 
integrand  in  the  integral  portion  of  the  cost  function  must  always  be  a  vector  quantity, 
both  for  the  Lagrange  and  the  Bolza  cost  functions. 

After  all  values  have  been  entered  the  file  is  ready  to  be  created.  When  the 


Figure  13.  Cart  Problem  Cost  Function  -  GUI 


Create  File  button  is  pushed  the  input  file  is  created  and  a  status  window  indicates 
the  task  completion,  see  Figure  14. 
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B.  OPTIMIZATION 

Once  the  file  has  been  created,  the  optimization  window  is  used  to  enter  the 
optimization  parameters.  The  name  of  the  file  to  be  optimized  is  entered  along  with 
the  number  of  LGL  points  to  be  used  in  the  optimization,  see  Figure  15.  The  larger 
the  number  of  LGL  points  the  finer  the  discretization.  Once  again,  fill  in  the  blanks 
before  pushing  any  buttons.  The  intuitive  flow  for  the  GUI  is  from  top  to  bottom, 
and  the  first  two  tasks  completed  need  to  be  the  selection  of  the  File  Function  Name 
and  the  Number  of  LGL  points.  The  number  of  LGL  points  may  be  changed  at  any 
time.  However,  changing  the  name  of  the  function  will  load  the  MATLAB  data  file 
associated  with  the  new  function  and  overwrite  the  current  variable  values. 


Figure  15.  Cart  Problem  Optimization  -  GUI 
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After  pressing  the  Initial  Guess  button  the  initial  guess  for  the  states  and  the 
controls  can  be  entered.  The  cart  problem  is  very  robust,  so  we  use  a  random  number 
generating  function  to  provide  a  random  initial  guess.  A  value  for  final  time,  t/,  is 
entered  along  with  the  initial  guesses.  If  tj  is  free,  this  is  an  initial  guess.  If  tj  is 
fixed,  the  value  for  tj  is  entered.  In  this  case,  the  fixed  final  time  value  is  tj  =  2,  see 
Figure  16. 


4  Initial  Guess 


•States, and  Cootrots:  A.  IniSaj  Guess” 

~  j'Y>-  ggi  m  -  *  ^ 


rand(1,n) 


Figure  16.  Cart  Problem  Initial  Guess  -  GUI 

The  options  for  the  constr.m  function  are  changed  by  the  Options  button. 
This  allows  us  to  set  upper  and  lower  bounds,  number  of  iterations,  and  number  of 
output  parameters.  We  can  also  check  the  number  of  equality  constraints  generated 
for  the  problem,  see  Figure  17.  There  are  four  groups  of  output  parameters.  The 
listbox  allows  the  selection  of  the  desired  set.  The  most  common  output  selection  is 
only  Xopt  because  the  addition  of  output  parameters  slows  the  optimization.  Xopt  is 
the  vector  of  optimization  variables:  the  states,  controls,  and  tj  for  a  free  final  time 


53 


problem.  The  value  of  the  MATLAB  output  parameter  vector,  the  problem’s  Hessian, 
and  lambdas  are  the  additional  output  parameters.  If  the  lambdas  were  the  Lagrange 
multipliers  for  the  NLP,  they  would  be  very  useful  for  costate  estimation.  Unfortu¬ 
nately  the  lambdas  returned  by  constr  .m  after  the  optimization  are  not  the  Lagrange 
multipliers  necessary  for  costate  estimation.  The  display  intermediate  results  check 
box  toggles  the  display  of  the  intermediate  steps  during  the  optimization  process. 

The  upper  and  lower  bound  can  be  left  blank  and  are  then  passed  to  constr. m  as 
empty. 


Figure  17.  Cart  Problem  Options  -  GUI 

Lastly,  we  input  the  values  for  the  parameters  designated  during  the  file  cre¬ 
ation,  see  Figure  18.  The  Direction  Vector  option  allows  designation  of  the  initial 
search  vector  direction.  The  default  is  set  to  empty.  The  Final  Time  figure  shows  a 
set  of  radio  buttons  to  check  the  status  of  the  final  time,  Figure  19.  If  this  button 
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Figure  18.  Cart  Problem  Parameter  Values  -  GUI 


Figure  19.  Cart  Problem  Final  Time  -  GUI 


does  not  agree  with  the  desired  value,  the  function  file  should  be  recreated  to  ensure 
the  proper  restriction  on  final  time. 

Once  the  optimize  button  has  been  pushed,  the  values  of  the  states  and  con¬ 
trols  will  be  output  to  the  MATLAB  command  window,  along  with  the  time  elapsed 
to  run  the  optimization,  and  the  final  value  of  the  cost  function.  In  addition,  an 
output  file  named  to  agree  with  the  input  file  is  created.  In  this  example  the 
output  file  is  cart  .function .  out.  Four  files  are  created  during  the  optimization: 
cart.f unction. m,  cart.functionopt  .a,  cart  .function. mat  and  cart.f  unction,  out. 

®  Is  the  m-file  that  contains  the  constraints,  state  equations  and  the 
cost  function.  The  cart .functionopt  .mis  a  script  that  will  declare  variables  and  call 
the  NLP  solver.  The  cart.f  unction,  mat  is  the  binary  file  that  saves  the  values  of 
the  variables  used  in  the  optimization.  The  cart.f  unction,  out  file  is  the  output  file 
that  captures  the  diary  of  the  output  of  the  optimization  to  the  MATLAB  workspace. 

The  GUI  implementation’s  results  agree  with  the  published  results  [Ref.  8].  The 
cart.f  unction,  m  and  cart  .function,  out  are  included  in  Appendix  A. 
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VII.  THE  ORBIT  TRANSFER  PROBLEM 


This  maximum  radius  orbit  transfer  problem  is  taken  from  the  text  by  Bryson 
and  Ho  [Ref.  14].  Given  a  constant-thrust  rocket  engine,  T  =  thrust,  operating  for 
a  given  length  of  time,  tj:  find  the  thrust  direction  history,  $(2),  to  transfer  a  rocket 
vehicle  from  a  given  initial  circular  orbit  to  the  largest  possible  circular  orbit.  The 
variables  are  defined  as 


r  —  radial  distance  of  spacecraft  from  attracting  center 


u  =  radial  component  of  velocity 

v  —  transverse  component  of  velocity 

m0  =  mass  of  spacecraft 

m  =  fuel  consumption  rate  (constant) 

$  =  thrust  direction  angle 

H  =  gravitational  constant  of  attracting  center 

The  problem  can  be  stated  as  finding  $(t)  to  maximize  r(tj )  subject  to  the  state 
equations 


v2  a  Tsm$ 

u  ~ - 2  - FTr 

r  rl  mo  —  \m\t 

uv  Tcos# 

v  = - 1 - FT# 

r  mo  —  |m|r 

with  initial  conditions 

r(0)  =  r0 
ti(0)  =  0 
®(0)  = 


(7.1) 

(7.2) 

(7.3) 
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and  final  time  conditions 


u(t/)  =  0 


This  problem  can  be  solved  with  fixed  final  time  as  a  maximum  orbit  distance  prob¬ 
lem,  or  with  free  final  time  to  find  the  minimum  time  to  a  specified  final  orbit.  Solving 
the  problem  as  a  maximum  radius  orbit  transfer  problem  with  fixed  final  time,  the 
cost  function  in  Mayer  form  is 

J  =  -r(tf).  (7.4) 

Equation  7.4  represents  maximizing  the  radius  by  minimizing  the  negative  of  the 
radius.  Before  finding  the  optimal  solution,  we  must  set  values  for  the  parameters  or 
constants  in  the  problem.  Let  • 

m0  =  1 
rn  =  -0.0749 
T  =  0.1405 
y  =  1 
r0  =  1. 

A.  FILE  CREATION 

The  implementation  for  this  problem  mirrors  the  cart  problem  of  Chapter  VI. 
The  state  equations  are  more  complex  than  those  required  for  the  cart  problem. 

r  =  u 

u  =  (v.~2)./r  -(mu) ./(r.~2)  +T*sin(phi) ./(m0-ml*t) 
v  =  (-u.*v./r  +  T*cos(phi) . /(mO-ml*t) 

are  entered  as  the  state  equations  into  the  State  Equations  GUI,  see  Figure  20.  All 
the  variables  in  the  state  equations  will  be  defined  in  the  creation  of  the  GUI,  with 
one  exception.  The  variable  t  is  predefined  within  the  GUI  to  represent  time.  The 
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Figure  20.  Maximum  Radius  State  Equations  -  GUI 

values  of  t  are  the  LGL  points  scaled  to  the  interval  [to,  t/].  After  entering  the 
controls  and  the  initial  conditions,  see  Figures  21  and  22,  the  final  time  conditions 
are  entered.  Figure  23  shows  the  final  time  boundary  conditions  when  final  time 
is  fixed.  This  problem  also  has  inequality  constraints  which  must  be  entered.  The 
inequality  constraints  are  restrictions  on  the  allowed  valued  for  the  control,  <j>.  The 
values  for  4>  are  limited  to  —7 r  <  4>  <  n.  Since  constr.m  only  allow  inequalities  of 
the  form  g(n)  <  0,  the  expression  must  be  transformed  into  the  normal  form  before 
it  may  be  entered  into  the  GUI. 
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Figure  21.  Maximum  Radius  Controls  -  GUI 


Figure  22.  Maximum  Radius  Initial  Conditions  -  GUI 


Figure  23.  Maximum  Radius  Final  Conditions  -  GUI 
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and 


—pi  <  <f)  <f>  —  pi  <  0 

pi  >  <j>  *->•  4>  —  pi  <0 


(7.5) 

(7.6) 


see  Figure  24.  The  constants  are  entered  in  the  parameters  GUI,  Figure  25.  The 


Figure  24.  Maximum  Radius  Inequality  Constraints  -  GUI 

final  step  is  the  creation  of  the  cost  function  as  shown  in  Figure  26.  The  discrete  cost 
function  for  the  Mayer  form  of  the  cost  function  is 

J  =  —  r(n). 

The  optimization  file  created  and  the  output  from  the  optimization  for  the  maximum 
radius  problem  are  located  in  Appendix  B. 
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Figure  26.  Maximum  Radius  Cost  Function  -  GUI 
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B.  OPTIMIZATION 

Before  optimizing  the  problem,  the  initial  guesses  must  be  entered,  Figure  27, 
The  initial  guess  for  each  of  the  states  and  controls  may  be  a  column  or  row  vector, 


Figure  27.  Maximum  Radius  Initial  Guess  -  GUI 

but  they  must  all  agree  in  dimension.  The  lin_intp  function  input  as  the  initial 
guess  for  r,  outputs  a  linear  guess  with  values  in  [1,1.525]  as  an  n  x  1  row  vector. 
Once  the  first  initial  guess  is  input  as  a  row  vector,  all  others  guesses  must  be  row 
vectors  also.  The  parameters  are  set  to  their  constant  values,  see  Figure  28.  The 
Options  button  can  be  used  to  check  the  number  of  equality  constraints  and  increase 
the  number  of  iterations  allowed.  The  Final  Time  button  can  be  used  to  check  for 
the  correct  time  characteristics,  Figure  29.  The  number  of  equality  constraints  in  the 
NLP  must  be  equal  to  the  number  of  state  equations  and  boundary  conditions. 
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Figure  28.  Maximum  Radius  Problem  Parameters  -  GUI 
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Figure  29.  Maximum  Radius  Problem  Options  -  GUI 
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Having  chosen  11  LGL  points,  the  number  of  equality  constraints  are 

3 n  =  33  state  equations 

3  initial/boundary  conditions 

2  final/boundary  conditions 

38  equality  constraints 

Often  the  default  number  of  function  evaluations  will  be  too  small  for  the 
problem  to  converge.  The  number  displayed  in  the  Options  GUI  is  the  default  for 
constr.m.  It  is  easily  increased  by  adding  zeros  to  the  displayed  number.  The 
optimize  button  is  pressed  to  begin  the  optimization.  A  diary  file,  bho.out,  resulting 
from  the  optimization  run  with  11  LGL  points  is  included  in  Appendix  B. 
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C.  CHANGING  THE  FILE 

Changing  the  problem  to  a  free  final  time  problem,  where  the  final  radius 
is  fixed  and  the  objective  function  switches  to  optimizing  time,  is  accomplished  by 
several  simple  changes  to  the  function  file.  By  pressing  the  Final  Conditions  button, 
adding  r(n)  -  radius  as  a  final  time  condition,  and  changing  the  final  time  to 
free,  the  problem  has  been  slightly  modified  to  the  new  variant,  see  Figure  30.  The 
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Figure  30.  Optimal  Time  Problem  Final  Conditions  -  GUI 

Inequality  Constraints  must  be  modified  to  include  a  non-negativity  constraint  on  tj, 
see  Figure  31. 

After  adding  the  parameter  radius  to  the  parameter  list,  Figure  32,  and 
changing  the  cost  function  to  tj ,  Figure  33,  the  new  problem  has  been  formulated. 
Pressing  the  Create  Function  button  creates  the  new  m-file  and  the  problem  is  ready 
to  be  optimized.  The  m-file  generated  is  included  in  Appendix  C. 
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Figure  31.  Optimal  Time  Problem  Inequality  Constraints  -  GUI 
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Figure  32.  Optimal  Time  Problem  Parameters  -  GUI 


Figure  33.  Optimal  Time  Cost  Function  -  GUI 


D.  A  SECOND  OPTIMIZATION 


Suppose  we  want  to  use  the  data  from  the  last  optimization  as  the  initial  guess 
for  solving  this  slightly  modified  problem.  Enter  the  values  from  the  last  optimization 
run  into  the  Initial  Guess  figure,  see  Figure  34.  Choose  a  value  for  the  parameter 
radius,  radius  =  1.345,  and  enter  the  value  in  Figure  35.  Check  the  Final  Time 


V  Initial  Guess  Efe  £3 


Figure  34.  Optimal  Time  Problem  Initial  Guess  -  GUI 

button  for  correctness  and  the  problem  is  ready  for  optimization.  Pressing  the  op¬ 
timize  button  optimizes  the  problem  and  appends  the  results  to  the  last  results,  see 
Appendix  D.  The  function  file  created  for  a  free  final  time  problem  no  longer  has  tj 
as  a  parameter  passed  into  the  function,  but  rather  assigns  the  value  during  run-time. 
Comparing  the  two  files  created  for  the  orbit  transfer  problem  illustrate  the  subtle 
differences  between  the  two  types,  fixed  and  free  final  time  problems. 
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Figure  35.  Optimal  Time  Problem  Parameters  Values  -  GUI 
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VIII. 


CONCLUSION 


The  MATLAB  GUI  for  the  Legendre-Gauss-Lobatto  Pseudospectral  algorithm 
makes  solving  optimal  control  problems  easier.  The  GUI  has  been  tailored  to  solve 
optimal  control  problems  in  MATLAB  using  the  Optimization  Toolkit.  This  provides 
a  wide  audience  of  potential  users.  Other  NLP  solvers  can  be  used  to  solve  trajectory 
optimization  problems,  such  as  NPSOL.  One  of  the  strengths  of  MATLAB  is  its  ability 
to  use  codes  written  in  non-MATLAB  languages  such  as  FORTRAN,  C,  or  C++. 
The  first  step  to  adapting  this  GUI  to  a  different  NLP  solver  is  to  understand  the 
inputs  that  the  solver  requires.  Once  the  input  parameters  have  been  determined,  the 
GUI  code  can  be  modified  to  create  different  inputs.  As  most  of  the  code  is  run  from 
functions,  the  GUI  figures  themselves  would  require  little  modification.  The  standard 
form  of  an  optimal  control  problem  will  not  change  with  the  addition  of  additional 
solvers.  The  Optimization  figure  would  require  the  most  visual  modification  to  add 
parameters  specific  to  each  solver.  The  GUI  could  be  easily  modified  to  create  several 
input  files  from  the  file  figure  and  then  allow  the  user  to  choose  from  several  solver 
options. 

The  GUI  can  be  used  to  allow  students  to  rapidly  solve  problems  in  this  area 
of  mathematics.  The  tool  can  also  be  used  to  provide  initial  guesses  for  numerical 
indirect  methods  that  require  very  good  initial  guesses.  The  GUI  has  been  beta 
tested  by  two  students  from  the  Aeronautics  and  Astronautics  Department,  LT  Bryan 
Schlotman,  USN,  and  Cpt  Lawrence  Halbach,  USAF.  Further  research  using  the  GUI 
will  be  conducted  by  these  two  officers  here  at  the  Naval  Postgraduate  School.  This 
GUI  will  also  be  used  in  the  Spacecraft  Performance  and  Optimization  (AA4850) 
class  in  the  fall  quarter. 

The  real  strength  of  the  LGLP  algorithm  involves  costate  estimation  from 
the  Lagrange  multipliers  returned  from  the  NLP  solver.  The  major  disappointment 
with  the  MATLAB  constr.mis  that  is  does  not  return  correct  Lagrange  multipliers. 
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Adaption  of  this  GUI  to  use  the  NPSOL  solver  will  allow  the  LGLP  method  to  be 
easily  used  and  evaluated  against  other  methods  using  the  de  facto  standard  sequential 
quadratic  programming  solver. 
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APPENDIX  A.  CART  PROBLEM 
OPTIMIZATION  FILES 


This  appendix  contains  one  of  the  files  created  while  solving  the  simple  cart 
problem  of  Chapter  6.  The  following  file  is  the  CONSTR.M  input  file  which  con¬ 
tains  the  function  to  be  minimized,  costfn,  and  the  constraint  function,  g.  The 
file  cart_f unction  .m  follows: 

function  [costfn, g]  =  cart_function(xopt,Dn,x,w,n,t,Tf  ,a,b,c)  ; 

^Define  the  variables 
xl=xopt(  0*n+l:n*  1)  ’ ; 
x2=xopt (  l*n+l :n*  2)’; 
u  =xopt(  2*n+l:n*  3)’; 


'/.Set  up  the  State  Constraints 
g(  0*n+l :n*  1)=(2/Tf )*(Dn*xl)-(x2  ); 

g(  l*n+l :n*  2)=(2/Tf )*(Dn*x2)-(-x2+u) ; 


%Set  up  the  initial  conditions 
g(  2*n+  l)=xl(l); 
g(  2*n+  2)=x2(2) ; 


%Set  up  the  final  time  conditions 
g(  2*n+  3) =a*xl (n)+b*x2 (n) -c ; 


V,Set  up  the  inequality  conditions 


’/.Set  up  the  cost  function 
costfn=(Tf/2)*(w,*(u.“2))  ; 
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APPENDIX  B.  MAXIMUM  RADIUS 
PROBLEM  OPTIMIZATION  FILES 


This  appendix  contains  two  of  the  files  created  while  solving  the  maximum 
radius  orbit  transfer  problem  of  Chapter  7.  The  first  file  is  the  CONSTR.M  input  file. 
The  file  bho  .m  follows: 

function  [costfn,gJ  =  bho(xopt ,Dn,x,w,n,t ,Tf ,mO,ml ,T,mu) ; 

y,Define  the  variables 
r  =xopt(  0*n+l :n*  1) ’ ; 
u  =xopt(  l*n+l :n*  2)’; 
v  =xopt(  2*n+l:n*  3)'; 
phi=xopt(  3*n+l:n*  4)'; 


#/,Set  up  the  State  Constraints 

g(  0*n+l:n*  1)=(2/Tf )*(Dn*r  )-(u  ) 

g(  l*n+l:n*  2)=(2/Tf )*(Dn*u  )-((v.~2)./r  -mu./(r.~2)  +  T*sin(phi) ./(mO-ml*t) ) 

g(  2*n+l:n*  3)=(2/Tf )*(Dn*v  )-(-(v.*u) ./r  +  T*cos(phi) ./(mO-ml*t)  ) 

%Set  up  the  initial  conditions 
g(  3*n+  l)=r(l)-l; 
g(  3*n+  2)=u(l)  ; 

g(  3*n+  3)=v(l)-l ; 

5lSet  up  the  final  time  conditions 
g(  3*n+  4)=u(n)  ; 

g(  3*n+  5)=v(n)-sqrt (l/r(n))  ; 

‘/.Set  up  the  point  inequality  constraints 

/.Set  up  the  inequality  constraints 
g(  3*n+  6:n*  4+  5)=phi-pi  ; 
g(  4*n+  6:n*  5+  5)=-phi-pi; 

'/.Set  up  the  cost  function 
costfn=-r(n) ; 
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The  bho .  out  file  resulting  from  an  optimization  run  follows: 


Filename_File  = 
bho 

date  = 

1999  6  7  17 

timer  = 

50.5900 

r  = 

Columns  1  through  7 

1.0000  1.0003  1.0033  1.0260  1.0898 

Columns  8  through  11 

1.3198  1.3397  1.3446  1.3450 

u  = 

Columns  1  through  7 

-0.0000  0.0046  0.0251  0.1034  0.1709 

Columns  8  through  11 

0.0821  0.0316  0.0085  -0.0000 

v  = 

Columns  1  through  7 

1.0000  1.0143  1.0418  1.0651  1.0028 

Columns  8  through  11 

0.8209  0.8320  0.8533  0.8622 

phi  = 

Columns  1  through  7 

0.6707  -0.0314  0.5439  0.4999  -3.1416 

Columns  8  through  11 
-1.1302  -0.5989  -0.7030  -0.6172 

cost  = 

-1.3450 


56  20 


1.1794  1.2660 


0.1871  0.1553 


0.8547  0.8158 


2.3196  -0.2570 
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APPENDIX  C.  OPTIMAL  TIME  PROBLEM 
OPTIMIZATION  FILES 


This  appendix  contains  two  of  the  files  created  while  solving  the  optimal  time 
orbit  transfer  problem  of  Chapter  7.  The  first  file  is  the  CONSTR .  M  input  file.  The  file 
bho .  m  follows: 


function  [costfn,g]  =  bho(xopt ,Dn,x,w,n,mO , ml ,T, mu, radius) ; 

Tf  =  xopt (length (xopt) ) ; 
t  =  (Tf/2)*(x+l) ; 

^Define  the  variables 
r  =xopt(  0*n+l:n*  1) ’ ; 
u  =xopt(  l*n+l:n*  2)’; 
v  =xopt(  2*n+l:n*  3)’; 
phi=xopt(  3*n+i:n*  4)’; 

%Set  up  the  State  Constraints 

g(  0*n+l :n*  1)=(2/Tf )* (Dn*r  )-(u  ) 

g(  l*n+l:n*  2)=(2/Tf )*(Dn*u  )-((v.~2)./r  -mu./(r.~2)  +  T*sin(phi) ./(mO-ml*t)) 

g(  2*n+l:n*  3)=(2/Tf )*(Dn*v  )-(-(v.*u) ./r  +  T*cos(phi) ./(mO-ml*t)  ) 

'/.Set  up  the  initial  conditions 

g(  3*n+  l)=r(l)-l; 

g(  3*n+  2)=u(l)  ; 

g(  3*n+  3)=v(l)-l; 

y.Set  up  the  final  time  conditions 

g(  3*n+  4)=u(n)  ; 

g(  3*n+  5)=v(n)-sqrt(l/r(n))  ; 

g(  3*n+  6)=r(n) -radius  ; 

/.Set  up  the  point  inequality  constraints 
g(  3*n+  7)— Tf; 

'/.Set  up  the  inequality  constraints 
g(  3*n+  8:n*  4+  7)=phi-pi  ; 

g(  4*n+  8:n*  5+  7)=-phi-pi; 
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y,Set  up  the  cost  function 
costfn=Tf ; 


The  bho.out  file  resulting  from  an  optimization  run  follows: 


Filename_File  = 
bho 

date  = 


1999 

6 

18 

13 

58 

27 

timer  = 

30.6500 

r  = 

Columns  1 

through  7 

1.0000 

0.9998 

1.0035 

1 . 0245 

1.0835 

1.1868 

1.2781 

Columns  8 

through  11 

1.3237 

1.3408 

1.3447 

1.3450 

u  = 

Columns  1 

through  7 

0.0000 

-0.0003 

0.0335 

0.0899 

0.1884 

0 . 2294 

0 . 1488 

Columns  8 

through  11 

0.0708 

V  = 

0.0288 

0.0059 

0.0000 

V 

Columns  1 

through  7 

1.0000 

1.0101 

1 . 0346 

1 . 0540 

1.0153 

0 . 8723 

0.7967 

Columns  8 

through  11 

0.8102 

0.8318 

0.8520 

0 . 8623 

phi  = 

Columns  1 

through  7 

-1.2427 

0.5160 

0.6096 

0.5087 

1.9920 

3.1416 

-0.9890 

Columns  8 

through  11 

-0.5805 

-0.6466 

-0.5589 

-0.3357 

cost  = 

3.1685 
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APPENDIX  D.  UTILITY  FUNCTIONS 


This  appendix  contains  the  utility  function  declare  .m.  This  file  is  used  within 
the  LGLP  GUI. 

function  [valstring, count]  =  declare (vars, vals) ; 
y.DECLARE  creates  variables  with  the  values  provided  in  the 
XMATLAB  workspace  (HI  line) 

'/. 

Xdeclare.m  -  A  utility  function  for  the  LGLP  GUI 

'/. 

7,  [valstring,  count]  =  declare(vars,vals) 

7. 

'/.input  vars  -  A  vector  of  strings  containing  the  names  of  variables 

'/. 

7%  vals  -  A  vector  of  strings  containing  the  values  to  be  assigned 

7,  to  the  input  variables 

'/. 

'/.output  valstring  -  A  vector  of  strings  that  when  executed  declare  variables 
7,  in  the  the  MATLAB  workspace.  The  command  ’ eval (valstring)  * 

'/.  will  execute  each  string  in  the  vector. 

'/.  count  -  A  counter  for  the  number  of  strings  or  size  of  valstring 

'/. 

'/. 

^written  by  Andrew  0.  Hall,  CPT,  US  Army,  June  1999,  at  the 
V,Naval  Postgraduate  School 

equal  =  ’  =  ’ ; 
semi  =  ‘ i ’ ; 
count  =  1 ; 
valstring  =  []  ; 
for  index  =  varsJ 

vstring  =  strcat (vars (count , : ) , equal)  ; 
vstring  =  strcat (vstring, vals (count , : )) ; 
vstring  =  strcat (vstring, semi)  ; 
valstring  =  strvcat (valstring, vstring) ; 
count  =  count  +  1; 
end 

count  =  count  -  1 ; 
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APPENDIX  E.  CONTENTS  FILE 


This  appendix  contains  the  Contents  .m  file  for  the  LGLP  GUI.  A  file  like  this 
needs  to  be  placed  within  a  directory  to  allow  the  file  to  be  accessed  through  the 
MATLAB  help  menus. 


X  LGLP  Toolkit  for  Optimal  Control  Problems 

'/.  written  by  Andrew  0.  Hall,  Naval  Postgraduate  School 

'/,  version  1.0  June  1999 

y. 

'/,  This  directory  contains  a  GUI  for  solving  optimal  control  problems. 

*/,  There  are  three  files  that  execute  the  GUI.  Typing  'opt’  at  the  command 
%  line  will  start  the  GUI. 

X 


X  GUI  Files 

X 

X  OPT.M 

X 


X  FILE.M 

X 

X 


X  OPTIMIZE. M 

X 


-  Welcome  screen  for  the  GUI  interface 

-  File  creation  GUI.  Creates  an  NLP  formulation  for  input  to 
CONSTR.M  (see  OPTIMIZATION  TOOLKIT) 

-  File  to  optimize  an  optimal  control  problem.  This 
function  sets  all  input  parameters  for  CONSTR.M. 


X  Utility  Files 

X 

X  costcall.m  - 

X 

X 

'/,  declare. m 

X 

X 

'/,  optcall.m 

X 

X 


y.  Example  Files 

X 

y,  carttest.m 

X 


Creates  a  function  call  to  compute  the  value  the  cost 
function 

Creates  a  function  call  to  create  variables  with  the  values 
provided  in  the  MATLAB  workspace 

Creates  a  function  call  for  CONSTR.M 


-  see  AOExamples.ps  in  this  directory 

-  Simple  caxt  problem,  linear  constraints  with  quadratic 
controls,  from  Chapter  6  of  Hall's  Thesis 
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V,  bh.m 
•/. 

V,  bho .  m 

y 


Maximum  Radius  Orbit  Transfer  problem  from  Chapter  7 
of  Hall's  Thesis 

Optimal  Time  Orbit  Transfer  problem  from  Chapter  7 
of  Hall's  Thesis 


%  Faculty  Codes 

y. 

y. 

y,  diffmat  .m 
*/.  lin.intp.m 
%  lobatto.m 

y. 

y,  mxt .  m 
y  mxtj  .m 
'/,  sgn.m 

y 

y 

*/.  tqr.m 


Codes  written  by  Professors  Bill  Gragg  and  Fariba  Fahroo, 
Naval  Postgraduate  School 

creates  a  differentiation  matrix  for  the  LGLP  algorithm 

creates  a  linear  guess  vector 

Computes  abscissa  and  weights  for  the  n-point 

Gauss- Jacobi-Lobatto  quadrature  formula 

creates  a  tridiagonal  matrix 

creates  a  Jacobi  matrix 

For  z  a  complex  number  we  define  sgn  z,  the  SIGNUM  of  z, 
as  z/|z|  if  z  =0  and  +  1  if  z  =  0.  Thus  sgn  z  is  the 
same  as  matlab's  sign  z  except  when  z  =  0. 

Tridiagonal  QR  algorithm 


y  Postscript  Files 

y 

y  AOExamples.ps  -  PS  file  that  includes  Chapters  5, 6, and  7  from  Hall's  Thesis 
*/.  A0Thesis.ps  -  Hall's  Naval  Postgraduate  School  Thesis  (Applied  Mathematics^ 


*/.  written  by  Andrew  0.  Hall,  CPT,  US  Army,  June  1999,  at  the 
'/.  Naval  Postgraduate  School 
y  Questions?  Email:  AndrewOHall@msn.com 
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